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Hit.  HI*  !•  (m  at  ih  1913;  **lHi  4 

dissertation  on,  "Asymptotic  Forms  of  the  Solutions  of  the 
Differential  Equation  for  the  Associated  Mathieu  Functions." 

He  began  his  college  teaching  career  in  1942  at  St.  Olaf 
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series  of  Forest  Products  Laboratory  reports  on  the  behavior  of 
sandwich  panels  under  loads,  and  continues  with  recent  articles  in 
the  Society  for  Industrial  and  Applied  Mathematics'  Journal  on  Nu¬ 
merical  Analysis,  on  inverse  pairs  of  matrices  with  integer  ele¬ 
ments. 
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THE  PCT  CONTROL  SYSTEM  DESIGN  FOR  SAMPLED-DATA  CONTROL  SYSTEMS 


C.  H.  Houpis* 


October  4,  1982 


ABSTRACT 


The  pseudo-conti nuous-ti me  (PCT)  control  system  approximation  of  a  sampled- 
data  system  permits  the  use  of  tried  and  proven  continuous-time  domain  methods 
for  designing  cascade  and/or  feedback  controllers.  When  the  rules  governing 
the  use  of  the  Pade  approximation  and  the  Tustin  transformation  are  satisfied, 
the  PCT  design  approach  is  a  valuable  technique  for  the  design  of  sampled-data 
control  systems. 


■'Professor  of  Electrical  Engineering,  School  of  Engineering, 

Air  Force  Institute  of  Technology,  Wri ght-Patterson  AFB ,  Ohio  45433 
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I.  Introduction 


The  analysis  and  design  of  sampled-data  control  systems  may  be  done  entirely 

in  the  z-plane,  which  is  referred  to  as  the  direct  digital  control  design  (DIR) 
f2l 

technique1-  J  or  entirely  in  the  s-plane.  The  latter  is  referred  to  as  the 
digitization  (DIG)  technique  which  requires  the  development  of  the  pseudo- 

I 

•'ontinuous-time  (PCT)  system  model.  This  model  requires  the  use  of  the  Pade 
and  Tusti  n  transformation  approximations.  A  control  ler^  0  (z)  designed  by  the  DIG 
technique  provides  a  good  base  for  exhibiting  the  effects  of  the  sample  time 
parameter  of  the  digitized  controller  on  the  performance  of  the  system.  The 
reason  for  this  is  that  the  continuous  controller  corresponds  to  the  limiting 
case  where  the  sampling  time  of  Dc(z)  is  zero.  A  disadvantage  of  this  method 
is  that  D  (z)  may  not  have  all  the  properties  of  D  (s).  However,  this  problem 

L  v 

is  minimized  by  the  selection  of  a  s  to  z  (to  w)  transformation  algorithm  which 
maintains  the  specified  properties  required  of  the  controller.  This  paper 
develops  the  PCT  control  system  model  and  the  criteria  for  achieving  a  good 
correlation  between  the  s-  and  z-plane  mapping.  If  the  degree  of  correlation 
is  not  satisfactory^  then  the  DIG  technique  may  permit  the  desired  system 
performance  characteristics  to  be  achieved  by  mere  gain  adjustment  in  the 


z-domai n. 


II.  Aoproximations 


Tustin  --  The  Tustin  s  to  z  or  z  to  s  transformation  is  defined  as 

saT<FT> 


,  .  T+sT/2 

z  *  ray? 


and  the  exact  z-transformation  is  defined  as 


By  substituting  s  =  Jw$p  into  C2J  and  s  «  j<2sp  into  (3),  where  u>sp  is  an 
equivalent  s-plane  frequency,  then  equating  the  two  expressions  results  in 


tan(u>spT/2)  =  w$pT/2 
When  d»spT/2  <  0.3057  rad  (17°)  then 


-  OJ 

sp  “sp 


-cr  T 


In  a  similar  manner,  substituting  the  exponentional  series  for  z  =  e  sp  and 

s  *  a,  into  (2)  results  in 
sp 

(a  T)2 

1  +  V  +  1  +  ...  =  1  +  [aspT/(l-aspT/2)]  (6) 

If  1  »  |aspT/2 1  and  1  »  f  crspT/2 1 ,  then 

l-spl  5  1‘spl  ~  2/T  (7) 

With  (5)  and  (7}  satisfied,  the  shaded  area  in  Fig.  1  represents  the  allowable 
location  of  the  poles  and  zeros  in  the  s-plane  for  a  good  Tustin  transformation. 

i  » 

Pade  --  Using  the  first-order  Pade  approximation,  the  transfer  function  of 
the  zero-order-hold  (Z-O-H)  device  is  approximated, when  the  value  of  T  is  small 

4 

enough  ,  as  follows; 


i  -  ? 
s 


JTT~7  s  Gpa(s) 


III.  Pseudo  Continuous -Time  Control  System  (PCT) 

The  DIG  method  of  designing  a  sampled-data  system,  in  the  complex  frequency 
s-plane,  requires  a  satisfactory  PCT  model  of  the  sampled-data  system.  In  other 
words,  for  the  sampled-data  system  of  Fig.  2,  the  sampler  and  the  Z-O-H  units 
are  approximated  by  a  linear  continuous -time  unit,  G^(s),  as  shown  in  Fig.  3(c). 
The  DIG  method  requires  that  the  dominant  poles  and  zeros  of  the  PCT  model 
must  lie  in  the  shaded  area  of  Fig.  1  for  a  high  level  of  correlation  with 


the  sampled-data  system.  To  determine  G^(s)  first  note  that  the  frequency 
component  of  E*(joi).  representing  the  continuous-time  signal  E(jw)  and  all  of 
its  side-bands,  is  multiplied  by  1/T  Because  of  the  low-pass  filtering 

characteristics  of  a  sampled-data  system,  only  the  primary  component  needs  to 
be  considered  in  the  analysis  of  the  system.  Therefore,  the  PCT  approximation 

l 

cf  the  sampler  and  the  Z-O-H  of  Fig.  3(a)  is  shown  in  Fig.  3(b)  where  the  Pade 
approximation,  Gpfl(s),  is  used  to  replace  Gzo ( s ) .  Therefore,  the  sampler  and 
Z-O-H  units  of  a  sampled-data  system  are  approximated  in  the  PCT  system  of  Fig. 
3(c)  by  the  transfer  function 


GA(s)  =  T  GpaCs)  =  TsVl 


Since  Lim  [G-(s)]  *  1,  (9)  is  an  accurate  PCT  representation  of  the  sampler  and 
T-fl  A 

Z-O-H  units,  satisfying  the  requirement  that  as  T  *  0  the  output  of  Gft(s)  must 
equal  its  input.  Further  note  that  in  the  frequency  domain,  as  ws  ■*■<*>  (T  -*■  0), 
the  primary  strip  in  the  s-plane  becomes  the  entire  frequency  spectrum  domain 
which  is  the  representation  for  the  continuous -time  system. 

Note  that  in  obtaininq  PCT  systems  for  the  sampled-data  systems  of  Fiq.  4. 


the  factor  1/T  replaces  only  the  sampler  that  is  sampling  the  continuous-time 


signal.  The  sampler  on  the  output  of  the  digital  controller  is  replaced  by  a 


factor  of  one.  To  illustrate  the  effect  of  the  value  of  T  on  the  validity  of 
the  results  obtained  by  the  DIG  method,  consider  the  sampled-data  closed-loop 
control  system  of  Fig.  2  where 

Gx(s)  3  sTs+Tffs+5T 

The  closed-loop  system  performance  for  three  values  of  T  and  c  =  0.45  are 
determined  in  both  the  s-  and  z-  domains,  i.e.,  the  OIG  and  DIR  methods, 
respectively.  Table  1  presents  the  required  value  of  and  time  response 
character!  sti  cs  for  each  value  of  T.  Mote  that  for  T  «-  0.1  there  is  a  high 


correlation  between  the  DIG  and  DIR  models.  For  T  <_  1  there  is  still  a  re¬ 
latively  good  correlation.  (The  designer  needs  to  specify,  for  a  given 
application,  what  is  considered  to  be  "good  correlation.")  The  figures  of 


merit  of  the  corresponding  continuous-time  control  system 
CM  .  Sx(s) 


for  a  unit-step  forcing  function  are:  =  1.202,  t  =  4.12  sec  and  t$ 

Table  1  Performance  Characteristics  of  a  Sampled-data  Control 
System  using  the  DIR  and  DIG  Methods 


(10) 

9.48  sec. 


Method 

T.sec 

«x 

mp 

Vsec  .  _ 

ts,sec 

DIR 

0.01 

4.147 

1.202  - 

4.16 

9.53 

DIG 

4.215 

1.206 

4.11 

9.478 

DIR 

0.1 

3.892 

1.202 

4.25 

9.8+ 

DIG 

3.906 

1 .203 

4.33" 

9.90+ 

DIR 

1 

2.4393 

1.2 

6 

13* 

DIG 

2.496 

1.200 

6.18 

13.76 

IV.  DIG  Technique 

Figure  5  represents  the  trial  and  error  design  philosophy  in  applying  the 
DIG  technique.  If  path  A  does  not  result  in  the  specifications  being  met  by 
the  sampled-data  control  system  of  Fig.  4,  then  path  B  is  used  to  try  to 
determine  a  satisfactory  value  of  KZ(;.  A  similar  chart  may  be  drawn  for  the 
design  of  the  feedback  controller  of  Fig.  4(b).  The  design  philosophy  involves 
the  following  considerations: 

(a)  Follow  path  A  if  the  dominant  poles  and  zeros  of  C( * )/R( - )  lie  in  the 
shaded  area  of  Fig.  1  (Tustin  aoproximation  is  good!). 

(b)  Follow  path  A  when  the  degree  of  warping  is  deemed  not  to  negatively 
affect  the  achievement  of  the  desired  design  results.  If  the  desired  results 
are  not  achieved^ try  path  B. 
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(c)  Follow  path  B  when  severe  warptng  exists.  The  DIG  design  procedure 
is  as  follows: 

Step  1  Convert  the  basic  sampled-data  control  system  to  a  PCT  control 
system  or  transform  the  basic  system  into  the  w-plane  (use  both  for  design 
Durposes) . 

Step  2  By  means  of  a  root-locus  analysis  or  by  use  of  the  Guilleman- 

I  i 

Truxal  method  determine  D  (sj  *  K  D  ($}  or  D  (w)  =  K  D  (w) . 

v  SC  L  L  Wj  C 

Step  3  Obtain  the  control  ratio  of  the  compensated  system  and  the  corre¬ 
sponding  time  response  for  the  desired  forcing  function.  (This  step  is  not 
necessary  if  the  exact  Guil lemin-Truxal  compensator  is  used.)  If  the  desired 
performance  results  are  not  achieved,  repeat  Step  2  by  selecting  a  different 
value  of  c.  <r,  etc.  or  a  different  desired  control  ratio. 

Step  4  When  an  acceptable  D  ("sj  or  D  (w)  has  been  achieved,  transform 
the  conpensator,  via  the  Tusttn  transformation,  into  the  z-domain. 

Step  5  Obtain  the  z-domatn  control  ratio  of  the  compensated  system  and 
the  corresponding  time  response  for  the  desired  forcing  function.  If  the 
desired  performance  results  for  the  sampled-data  control  system  are 
achieved  via  path  A  or  path  B,  then  the  design  of  the  compensator  is  complete. 
If  not,  return  to  Step  2  and  repeat  the  steps  with  a  new  compensator  design 
or  proceed  to  the  DIR  technique. 

V.  DIR  Technique 

The  simple  lead  fa  <  1 )  and  lag  fa  >  1)  compensators  in  the  s-domain  have 
the  form 


0c(s)  ■ 


Kscfs  ' 

(s  -  Pj 


where  ps  3  z^/i.  The  s-plane  zero  and  pole  are  transformed  into  the  z-domain 
as  follows: 


AAUl'A'. 
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■  ■  ■  J 


zz  =  f' 


pz  -  e- 


PST  ZsT/a 

=  f  -  e 


Thus,  the  corresponding  first-order  z-domain  compensator  (digital  filter)  is 

n  .  K*C(Z  -  .  "zc'2  -  2Z> 

V2'  rz  z/s 

where  pz  =  z z/s. 

By  taking  the  natural  log  of  (12)  and  (13)  a  relationship  between  a  and  0  is 
obtained  as  follows: 


V  =  1n  zz 
z$T/a  -  In  pz 


Taking  the  ratio  of  these  equations  and  rearranging  yields  a  In  n  =  In  z  . 


"z  ■  2z 


8  *  Zz/P2  ■  P, 


For  a  lead  network,  a  <  1  and  pz  is  also  less  than  one.  Therefore,  for  a 
lead  digital  filter  3  >  1.  For  a  lag  digital  filter  e  <  1.  (Note  that  the 


condition  on  3  is  just  the  opposite  for  that  on  a.) 

Because  (11)  and  (14)  have  the  same  mathematical  form , 


the  z-plane 


compensator  design  procedures  via  the  DIR  techniques  are  essentially  the  same 
as  those  for  designing  a  compensator  for  a  continuous-time  system. ^ 


VI.  Example  of  DIG  Design:  Guil laman-Truxal  (G-T)  Compensation  Method 


The  figures  of  merit  for  the  control  system  of  Fig.  2  where 

r  /.i  .  0.4767 
Gx(s)  ‘  sTi+TJ 

and  T  =  0.1  sec  are:  Mp  =  1.043,  t  =  6.45  sec,  t$  =  8.65  sec,  and  K]  =  0.4765  sec' 
The  control  ratio  for  the  corresponding  PCT  system  is 


CCs)  , 

mr 


9.534 


W.' 


Consider  the  case  where  the  figures  of  merit  of  the  basic  system  are  to  be  improved 
as  follows:  tp  and  t$  are  to  be  cut  by  one-half  and  some  improvement  in  K-j  is 
desired  while  maintaining  <_  1.10.  Further,  assume  that  the  compensator  model 
of  Fig.  4(a)  is  constrained  to  increase  the  order  of  the.  system  to  four. 

Based  upon  these  specifications  the  following  factors  are  used  to  derive  the 
desired  control  ratio  model: 


(1)  The  real  part  |<jj  2I  -^e  dominant  roots  is  selected  to  be  at  least 
twice  that  of  (19)  based  upon  T$  =  4/h,  2 | >  (2)  the  dominant  roots  are  selected  such 
that  c  =  0.7071  in  order  to  try  to  maintain  Mp  <_  1.10,  and  (3)  the  s-plane 
pole-zero  combination  of  z-j  =  -1.4  and  p^  =  -1.1  is  added  to  minimize  the 
increase  in  the  overshoot  which  occurs  when  transforming  from  the  conti nuous¬ 


time  model  to  the  sampled-data  model.  This  selection  of  values  is  made  in 


order  to  meet  the  desired  performance  specifications.  Thus^the  following 
continuous-time  control  ratio  model  is  achieved. 


15  714(s  +  1.4) _ 

(sZ  +  2s  +  2}(s  +  l.l)(s  +  10) 


(20) 


Although  the  zero-pole  combination,  -1.4  and  -1.1,  of  (20)  lie  just 
outside  the  allowable  region  of  Fig.  1,  this  aspect  is  overlooked  for  a  first 
trial  design.  Applying  the  Tustin  transformation  to  (20),  for  T  =0.1,  yields 


Applying  the  G-T  method  yields  the  following  transfer  function  of  the 


digital  compensator  of  Fig.  4. 


E,(z) 


of  Dc(z)  by  the  d-c  gain  factor  of  2. 

VII.  Conclusions 

By  the  use  of  the  Pade  approximation  and  the  Tustin  transformation,  a 
sampled-data  control  system  may  be  transformed  into  a  PCT  control  system.  As 
the  examples  illustrate,  when  the  rules  governing  this  transformation  are 
satisfied  the  analysis  and  design  of  a  PCT  system  model  is  a  practical  approach 
for  the  analysis  and  design  of  a  sampled-data  control  system.  The  standard 
first-order  z-plane  compensator,  Kzc(z-zz) / (z- zzfi),  corresponds  to  the  standard 
first-order  s-plane  compensator,  Ksc(s-zs)/(s-zs/a) ,  where  for  a  lead  compensator 
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i  <  I  and  a  >  1.  A  consequence  of  applying  the  Tustin  transformation  to 
C(s)/R(s)  of  the  PCT  system  and  then  applying  the  G-T  method  to  [C(z)/R(z)]TU 
is  that  an  unrealizable  cascade  digital  compensator  Dc ( z )  results.  This  paper 
illustrates  a  method  by  which  this  Dc(z)  may  be  made  realizable  by  replacing 
one  or  more  (z+1)  factor  in  the  numerator,  due  to  the  Tustin  transformation, 
by  its  d-c  gain  factor  of  two. 
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compensated  samp led- data  or  a  digital  control  system: 


0(0  -  M*)c»(z) 

R(z)  l  +  D c(z)CA(z) 


DIG  design  philosophy  for  Fig.  U-i(*) 
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ASYMPTOTIC  NON-NULL  DISTRIBUTION  OF  A  TEST 
OF  EQUALITY  OF  EXPONENTIAL  POPULATIONS 


R.  W.  Kulp  and  B.  N.  Nagarsenker 


Air  Force  Institute  of  Technology 
Wrlght-Patterson  Air  Force  Base,  Ohio 


Key  Words  and  Phrases:  exponential  populations;  likelihood  ratio  criterion; 
asymptotic  non-null  distribution;  Chi-square  distributions. 


Abstrac  t 


In  this  paper  asymptotic  expansions  of  the  non-null  distribution  of  the 
likelihood  ratio  criterion  for  testing  the  equality  of  several  one  parameter 
exponential  distributions  are  obtained  under  local  alternatives.  These  expan¬ 
sions  are  in  terms  of  Chi-square  distributions. 


Introduction 


Suppose  that  p  samples  are  available  and  that  the  ith  sample  contains  n 
observations  x  with  mean  x^  (i  *  l,2,...,p;  j  =  l,2,...,n)  and  has  been  drawn 


from  an  exponential  distribution  with  probability  density  given  by 


f(x)  •  oi  exp  (-x/o1> 


x  >  0,  oL  >  0 


otherwise  ( i=l , 2 , . . . ,p) 


(1.1) 


A  test  of  hypothesis  that  the  p  samples  have  been  randomly  drawn  from  the 
same  population  is  equivalent  to  testing  that  the  p  exponential  distributions 
in  (l.l)  are  identical.  In  other  words,  it  is  desired  to  test  the  hypothesis 


V  °1  =  °2 
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gainst  the  general  alternatives.  The  likelihood  ratio  criterion  for  testing 


Proof .  To  obtain  the  hth  moment  of  we  shall  essentially  use  the  method 
given  in  Wilks  (1946)  and  the  fact  that  are  independently  distributed  as 


gammas  with  parameters  n  and  ?Jn,  i  *  l,2,...,p.  For  this  consider  the 
function  where 

HO)  -  E  [  |  TT  <;.)nh  I  I  . 

*1*1  1  1  J 

It  can  be  easily  shown  that 

«0>  -  [  ]  PlE/n|nh|l-eiZr(nh+n)  (2.2) 

where  9^  *  9/np.  Using  the  following  identity  (see  Khatri  and  Srivastava  (1971)) 

|l-9  E|"(n+nh)  -  (l-91q)"p(nh+n)|q-1rr(nh+n)  £  £ 

'  1  ~  k*0  <  (l-9iq)k  k! 

where  0  <  q  <  °»  and  can  be  chosen  such  that  the  expansion  in  the  series  form  is 
valid,  i.e.,  9^  and  q  are  such  that 

^■q^bmaxP'1!  < 

E(lb)  is  then  obtained  by  evaluating  at  9  *  0  and  then  putting  r  *  -nh 

d9r 

(see  Wilks  (1946)  for  the  validity  of  such  operation).  This  gives  (2.1). 

2 

Remark.  Taking  q  *  o  ,  we  get  the  null  moments  of  X  given  in  Nagarsenker  (1980). 


SI 


Lemma  2 .  Let  C^(Z)  be  a  zonal  polynomial  corresponding  to  the  partition 

ic  ■  (k,  ,k„, . . .  ,k  }  with  k,  +  k_  +  ...  +  k  =  k  and  k,  >  k_  k  ...  >  k  >0.  Let 
12  p  12  p  1Z  p 


aL(K)  *  £  k.(kt-i),  a2(<)  *  [  k . (4k  2-6ik  +3i2) , 

i*l  1  i=l  1  1 


and  0^  *  tr(Zr).  Then  the  following  equalities  hold: 


jtotic  Non-Null  Distribution  of  L 


In  this  section  we  shall  obtain  the  asymptot ic  expansion  ?f  the  distributic 
of  -2m  In  L  in  terms  of  m  *  n  -  u  increasing  where  u  =  (p+l)/6p  (see  Nagarsenker 
(1980)),  for  the  sequence  of  alternatives  stated  in  (1.3). 
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Case  1.  (I-ql  ^ )  =  ?/m. 

Let  U  *  -2m  In  L.  Then  from  (2.1)  the  characteristic  function  Ht)  of  U  under 
this  sequence  of  alternatives  is  given  by 


Ht)  -  (  P-2»H  ]  P  l  l 


«  HtcK^P/m) 


k=0  k: 


k! 


(3.1) 


where  g  *  ( 1—2 it)  and  T(<) 
Now  using  the  expansion 


(mg+u) (pm+pu+k) 
r (pmg+pu+k) 


2  3  —4 

ln|l-P/m|  *  -o^/m  -  O^/^m  -  a^/3m  +  0(m  ) 


where  ■  tr(P  ),  we  have 


I-P/m | 


(m+u) 


~°1  f  Ci  C2  _3  l 

e  1  [  1  +  —  +  -j  +  0(m  J)  | 


(3.2) 


m  m 


where  C. 


- -  ucy  and 


Again  using  Stirling's  asymptotic  formula  for  the  logarithm  of  a  gamma  function 

_3 

and  then  using  lemmas  1-3  and  (3.2)  in  (3.1),  we  have  up  to  0(m  ), 


*(t)  *  g_v|  1  +  2^  (1-8_1)  +7T~72  l  3igl  +  0(m~3)l 

r  (2pm)  i“0 


(3.3) 


where  v  *  ,  b]_  =  ~  P°2*  b2 


ai  -  P  a3; 


SQ  *  y  b^2  +  y  b2  +  2uPb1  +  (p+l)2(p-l)/36 


Si  *  -b^  -  LbL  -  2b^  -  Aupbj^ 


and  S. 


“(w- 
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7?",  T7  w  4J"  *J~  *J~  ’  .  S  /«  ‘  „*V 


L  .•■•>. 


Inverting  the  characteristic  function  _ (t)  in  (3.3),  we  have  the  follouir 


theorem. 


Theorem  3.1.  Under  the  sequence  of  alternatives  I-ql'1  =  P/m,  the  non-null 
distribution  of  -2m  In  L  where  L  is  given  in  (1.2)  can  be  expanded  asymptotically 
for  large  m  -  n  -  u,  u  -  (p+l)/6p  as  follows: 

2  k.  y  j 

P(-2m  In  L  <  x)  »  P(X  <  x)  +  [ P (X  _  <  x)  -  P(X 2  ,  <  x)  ] 

t  2pm  f  f+2 


2 

+  ~  2  [  l  3iP(x2f+2i  -  x)|  +  0(m"3) 

(2pm)Z  i-0  1  £  21  1 


where  f  *  2v  3  p  -  1,  X  ^  is  the  chi-square  variable  with  f  degrees  of  freedom 
and  the  coefficients  b^ ,  8^,  8^  and  ^  are  given  in  (3.3). 


Case  2.  ( I— q  I)  *  Q/m. 

We  have  I-q£  ^  ■  -Q(I-Q/m)  £/m.  So  proceeding  as  in  the  case  1  by  replacing 
P  by  -Q(I-Q/m)  and  retaining  terms  of  the  order  of  m  we  have  the  following 
theorem. 


Theorem  3.2.  Under  the  sequence  of  alternatives  I  -  q  Z  =  Q/m,  the  non-null 
distribution  of  -2m  In  L  can  be  expanded  asymptotically  for  large  m  as  follows: 


P (-2m  In  L  <  x)  =  P(X2  <  x)  +  — -  [P(X2  <  x)  -  P(X2  -  i  x) 

f  2pm  f  f+2 


it?  I  ^  a‘p«2,«i  *- *»  I  +0(""3> 


2  2 

where  a^  *  (tr  Q)  -  p  tr  Q  , 


3  2  3 

a2  *  (tr  Q)  -  p  tr  Q  , 


■*  2>‘  / 1'.'*1  *' 


20  =  (p+1) 36P-1)  +  7  al2  "  4(tr  9)al  +  2uPai  +  3a?/3 


=  -a  +  2a^  +  8  (tr  Q)a^  -  4upa^  +  ia2 


and  a2  *  -(dg+a^) 


Remark.  It  may  be  noted  that  when  P  ■  Q  *  0,  the  asymptotic  expansion  in  the 


two  cases  reduces  to  that  of  Box’s  approximation  given  in  (4.2)  of  Nagarsenker 


(1980)  for  the  null  hypothesis  H^. 
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Irrotational  and  Solenoidal  Waves 
in  General  Coordinates 

by 

D.  A.  Lee 

Air  Force  Institute  of  Technology 
Abstract 

Tensor  analysis  and  some  identities  are  applied  to  make  tolerably  convenient 
forms  of  equations  for  irrotational  and  solenoidal  components  of  elastodynamic 
displacement  fields,  valid  in  arbitrary  admissible  curvilinear  coordinates. 


Intro  due  tion 

Displacement  fields  of  elastodynamic  waves,  i.e.  suitably  smooth  vector 
functions  u(x,t)  which  are  solutions  of  the  Navier-Cauchy  equations 


azV(VK)-6zVxV>m=!p£-  (l) 

can  always  be  written  as  sums  of  irrotational  and  solenoidal  vector  functions,  in 
the  form 


li=V^(y,z,f  )+7)t4 
where  the  vector  function  A(z..t)  is  solenoidal,  i.e. 

V.4=0 

(Reference  [  1  ]).  The  scalar  satisfies  the  scalar  wave  equation 

ozvV=^£- 
r  9t 2 

and  the  vector  A  satisfies  the  vector  wave  equation 


In  (5). 


V.*4=V(V  -A) -VxVxji 


(2) 

(3) 

(4) 

(5) 


in  view  of  (3)  (Reference  [2]). 


=  -V  xVxA 


Often  it  is  more  convenient  to  find  the  scalar  <p  and  the  vector  .A  represent¬ 
ing  an  elastodynamic  displacement  field,  than  it  is  to  find  directly  the  field  itself. 
Often,  too.  one  wishes  to  work  in  a  coordinate  system  convenient  for  the  prob¬ 
lem  at  hand,  and  then,  particularly  for  non-orthogonal  curvilinear  coordinates,  a 
blunt  evaluation  of  the  spatial  differential  operators  in  (4)  and  (5)  may  be  very 
unwieldy.  This  paper  presents  a  development  of  the  representation  (2)  valid  for 
arbitrary  curvilinear  coordinates,  together  with  convenient  means  for  writing 
out  generalizations  of  (4)  and  (5).  In  particular,  the  generalized  equations  can 
be  written  out  explicitly  without  evaluating  Christoffel  symbols. 


I rrotational  and  Solenoidal  waves  in  General  Coordinates 


In  general  coordinates,  an  irrotational  wave's  displacement  field  can  be 
written 

(9) 

while  a  solenoidal  wave's  displacement  field  can  be  expressed  in  the  form 

ui=-j==-ti>k 

where  the  once-covariant  tensor  4k  is  a  function  of  position  and  time,  as  is  the 
scalar  y  . 

In  (9)  and  (10),  g y  denotes  the  contravariant  metric  tensor  elements  of 
the  coordinate  system.  A  comma  preceeding  an  index  denotes  covariant 
differentiation,  e.g. 


(10) 


04k 


An 


0*T 

The  permutation  tensor  elements  may  be  defined  as 

[ spn  ( rr)  if  (i.j.fc )  is  a  permutation  n  of  (1,2.3) 
0.  otherwise 


=e 


u* 


These  elements  are  special  in  that  they  transform  both  as  three-times  contra¬ 
variant  relative  tensor  elements  of  weight  1  and  as  three-times  covariant 
relative  tensor  elements  of  weight  -1  . 

The  Navier-Cauchy  equations  (4)  may  be  written  in  general  coordinates  as 

o  zukJtm9mi-f>2eiikekrmumJfjgi,r=~~-  (11) 

Seeking  an  irrotational  solution  of  ( 1 1).  one  finds  that 

umP9pk=gmr?rP9’* 

is  symmetric  in  m  and  k,  so  that  the  second  term  of  the  left  side  of  (11)  is  ident¬ 
ically  zero.  The  remaining  terms  of  (11)  require 


which  will  surely  be  met  if 


(12) 

One  may  take  (12)  as  the  defining  equation  for  irrotational  waves. 

Turning  to  solenoidai  waves  with  displacement  fields  of  the  form  (10).  one 
sees  that 


<1/ m  = 

u  m 


0 


by  the  symmetry  of  Jm  with  respect  to  j  and  m,  so  that  now  the  first  term  on 
the  left  side  of  (ll)  is  identically  zero.  The  remaining  terms  lead,  after  some 
manipulation,  to 


a2^ 


af! 


=o 


(13) 


for  which  it  is  sufficient  that  the  quantity  in  square  brackets  is  identically  zero. 
But  that  condition  leads,  after  use  of  the  identity 


r 

n 


to 


But  the  defining  tensor  A*  of  the  solenoidai  waves  is  itself  solenoidai.  so  that 

g^A,j,=0  (14) 

Then  the  defining  equations  for  a  solenoidai  wave  may  be  written  as  (14), 
together  with 

(15) 

One  could,  of  course,  have  written  (12),  (14)  and  (15)  as  tensorially  con¬ 
sistent  generalizations  of  (4), (3).  and  (5),  respectively.  It  is  well,  however,  to 
carry  through  a  complete  treatment,  starting  with  explicit  definitions  (9)  and 
(10),  to  be  certain  the  work  is  self-consistent. 

Moreover,  while  the  forms  of  tensor  equations  (12)  and  (15)  are  succinct 
and  easy  to  remember,  these  forms  aren’t  usually  the  most  convenient  ones  to 
use  for  writing  out  the  detailed  statements  which  those  equations  imply  in  a 
given  coordinate  system.  Direct  evaluation  of  (15),  for  example,  usually 
requires  evaluation  of  ChristofTel  symbols  and  derivatives  of  ChristofTel  symbols. 

Explicit  ''unfolding''  of  (12)  for  a  specific  coordinate  system  is  made  easier 
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by  the  well-known  identity 


_JL _ a_ 

ax' 


by  virtue  of  which  (12)  is  equivalent  to 


d<fi  -dz» 


Vp  azn  (  *  a*n  J  a*8 

Explicit  writing -out  of  (15)  is  probably  simplified  in  most  cases  by  the  fol¬ 
lowing  considerations:  The  elements  Ti,  defined  by 

are  once-contravariant  (oriented)  tensor  coordinates.  Also,  they  may  be 
evaluated  explicitly  without  evaluating  ChristofTel  symbols,  since 

(l8) 

and  the  last  term  in  (18)  is  identically  zero  by  virtue  of  the  symmetry  of  the 
ChristofTel  symbol  in  its  lower  indices  j  and  k.  and  the  anti-symmetry  of  z1**  in 
those  indices.  Thus  the  Si. 

S  r  =  —  zr'k 

are  once  covariant,  oriented  tensor  coordinates.  Repeating  the  previous  argu¬ 
ments  of  this  paragraph  then  shows  that  the  Rit 

Spn 

—  Sim  cmnr>  _  Sim  gmnp  d  I  Spr  £ rjk  ^^k  1 


are  once  covariant  tensor  coordinates,  which  can  be  written  out  without  evaluat¬ 
ing  ChristofTel  symbols.  In  rectangular  Cartesian  coordinates. 


Rk  =Ekvp  Bpjr 


dx>dxn 


Now,  the  quantity 

s  £fcnm  Emr*  \.rp 

appearing  in  (13)  is  also  a  set  of  once  covariant  tensor  coordinates.  In  rectangu¬ 


lar  cartesian  coordinates, 


Kit  — 


d2\ 


But  then  =  ,  since  coordinates  of  a  given  tensor  character  which  are  equal 

in  one  coordinate  system  are  equal  in  all  coordinate  systems. 

Thus  a  set  of  equations  equivalent  to  (14),  (15)  is  (14)  and 


2  P  3frr  ^  ^<4^  _  £  A 

Vg  dxn  'Sg  dz>  dt2 


(19) 


Equation  (19)  is  likely  to  be  more  convenient  for  writing  out  than  equation  (15), 
since  no  Christoffel  symbols  appear  in  (19).  Moreover,  the  operations  on  the  left 
side  of  (19)  can  all  be  accomplished  by  fairly  straightforward  matrix  multiplica¬ 
tions  and  differentiations. 
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SOME  PROBLEMS  OF  ESTIMATION  WHEN  SOME 


PRIOR  INFORMATION  IS  AVAILABLE 

Albert  H.  Moore 

Air  Force  Institute  of  Technology 
Wright -Patterson  Air  Force  Base 

Abstract 

In  this  paper  we  will  consider  problems  where,  without 
prior  information,  the  classical  best  unbiased  estimators  are 
known  to  be  admissible.  If  we  assume  enough  prior  information 
is  available  to  bound  the  parameter,  but  not  enough  to  specify 
a  prior  density  we  will  show  that  the  classical  U.M.V.U.  estima¬ 
tors  are  inadmissible  by  exhibiting  alternative  biased  estimators 
with  uniformly  smaller  mean  square  error.  This  is  possible 
because  if  we  examine  the  mean  square  error  we  see  that  the 
uniformly  minimum  variance  unbiased  (U.M.V.U.)  estimators  are 
best  only  at  the  boundary  points  in  the  parameter  space.  If 
enough  prior  information  exists  to  exclude  these  boundary  points 
Chen  the  classical  estimator  should  be  inadmissible.  With  this 
insight  we  are  able  to  show  that  the  classical  U.M.V.U.  esti¬ 
mators  are  inadmissible  in  a  new  way. 

In  Sections  2-6,  we  give  several  examples  where  the  U.M.V.U. 
estimator  is  inadmissible  when  enough  prior  information  is  avail¬ 
able  to  bound  0.  Throughout  this  paper  we  use  squared  error 
loss  function. 


29 


1.  Introduction 


Let  x  be  a  random  variable  taking  values  in  a  measure 

space  (X,  B,  y)  .  In  standard  statistical  problems  X  is  an 

n-dimensional  Euclidean  space  E  and  B  is  a  Borel  a-field 

n 

and  |i  is  either  a  counting  measure  or  Lebesgue  measure.  It 
is  also  assumed  that  x  has  a  density  f(x,6),6eft  *'  parameter 
space,  with  respect  to  the  appropriate  measure. 

In  estimation  problems,  we  assume  f  is  known  but  6  is 
unknown.  We  observe  X  »  x,  and  want  to  find  a  measureable 
function  of  X  *  x,  and  use  it  to  estimate  6.  The  classical 
approaches  are  maximum  likelihood  estimation  (M.L.E.)  and 
li.M.V.U.  estimation.  In  the  decision  theory  approach  Bayes 
and  minimax  estimators  are  generally  used. 

In  U.M.V.U.  estimation  any  prior  information  is  disre¬ 
garded  while  on  the  other  hand  in  many  instances  the  Bayes 
approach  requires  too  much  prior  information.  In  many  cases 
we  have  only  partial  information. 

Robbins  [5]  proposed  an  empirical  Bayes  approach  which 
makes  use  of  past  data  to  obtain  an  empirical  estimate  of 
the  prior  density  of  the  parameter  for  Bayes  estimation. 

Katz  [3]  investigated  some  properties  of  point  estimators 
when  an  upper  or  lower  bound  for  the  parameters  is  given  in 
advance.  He  considered  the  square  error  loss  function  and 

derived  adml  *  s  t  M  I  ti  ■*  n  -  it  <  •i*  -  ••-‘•■p  "-.  *-<■*>«  -ft'*"1'*' 

with  a  suitable  prior.  Lehmann  [4]  showed  that  for  N  ( 9  , 1 ) 
with  a  <  9  <  bs x  is  not  admissible  and  not  minimax.  In 


addition  he  showed  that  if  9  >  a,  then  x  is  minimax  but  not 


admissible.  Skibinsky  and  Cote  [6]  showed  that  with  certain 


prior  information  about  the  distribution  of  9  for  the 


binomial  x  is  inadmissible,  and  in  a  similar  fashion  he 


showed  that  x  is  an  inadmissible  estimator  for  the  mean  of 


the  normal  density.  Kale  [2]  showed  that  for  truncated 


parameter  spaces  the  M.L.E.  is  inadmissible  by  showing  it 


is  not  a  proper  Bayes  estimator  for  the  exponential  family 


of  densities  with  continuous  parameter  space.  Blum  and 


Rosenblatt  [1]  discussed  Bayes  estimation  where  it  was 


assumed  that  the  family  of  distributions  which  the  prior 


came  from  is  known. 


2.  The  Binomial  Distribution.  Let  X  be  an  observa- 


-tion  from  a  binomial  distribution  with  density  b(x,n,8)  ■ 


*  0x(l-8)n  x,  0  S  8  i  1.  The  statistic  X  is  a  complete 


sufficient  statistic  for  9.  The  best  unbiased  estimator  6, 


of  0  is  x/n.  The  risk  of  8^  is  given  by 


(e)  .  Uiliill  . 

e,  n 


Consider  an  estimator  6^  of  the  form  6^  »  k9 ^ .  The 


value  of  k  which  minimizes  the  mean  square  error  is  easily 


seen  to  be 


.  q" 

n 


L*±  ♦  e 

r» 


sc 


Here  the  best  estimator  of  the  fora  kQ ,  is  a  function  of  the 
unknown  parameter  0.  However  we  can  say  that 


TX— 7  V°  s  ‘  5  1 

IT  1 


is  the  best  estimator  for  9  of  the  form  k8^  at  the  parameter 

value  0  ■  A.  It  is  interesting  to  see  that  the  U.M.V.U. 

estimator  is  best  (in  the  mean  square  sense)  only  at  the 

parameter  value  0  ■  1.  The  estimators  6^  are  an  infinite 

family  of  admissible  estimators.  Suppose  we  have  prior 

information  that  0  i  t  <  1. 

Theorem  1.  If  it  is  known  that  0  $  l  <  1  then 

*  *$i 

0j  ■  jTj -  is  an  estimator  with  uniformly  smaller  mean 

n  *  1 

^  A 

square  error  than  0^  and  hence  0^  is  inadmissible. 


Proof.  Consider  9,  ■  k9,  where  k  • 


(l-l)  ♦  ni  * 


......  of  $  .  ■.V(l-9).2nt9.nt29.n9 

“  (1-0)  (l-2JU2nt-2nt  -n2lZ] 


eq-9) 


n0ri-21*&2l-*-n2t2(l-8) _ 

(l-0)[l-2£*i2]-»'2nl(l-0)*n2i2(l~9) 


n0  ♦ 


9(1-9) 

n 


l-2i*l2 


(1-0)  [l  ♦  UUdlzU  ♦  Slihlli 

L  (1-t)  1-21*1* 


32 


e  -  m-a 1 


_ n9*n2l2i;i-9)/l-21<-l2 _ 

2(i-t)9)  *  (1"9:  4  n2t2(l-9)/l-2i*t2 


Since  9  4  4  the  fraction  2n^^~9-  >  n9 .  Therefore 


Re  <®>  <  6(n~^  *  Re  <9> 
03  n  81 


and  hence  9 ^  is  inadmissible. 


3  .  The  Poisson  Distribution.  Consider  the  Poisson 


density 


e‘AXx 

f(x,X)  -  *  ,  x  -  0,  1,  2,  ...  .  (4) 

Let  x..  ....  x  be  a  random  sample  from  a  Poisson 
1  "  n 

distribution.  Let  X.  ■  J  x./n.  X.  is  a  U.M.V.U.  estimator 

1  i-1  1 

of  X  with  risk  R j  (X)  *  X/n.  Consider  the  estimator  X.  of 

*1  * 

A  A  * 

the  form  X2  *  kXj.  The  value  of  k  which  minimizes  the  mean 
square  error  is  easily  seen  to  be 


X4  ♦  X/n 


Therefore  an  estimator  of  the  fora 


X.  -  - = -  X  .  I  >  0 

3  l  ♦  /it  n  1 


is  the  best  (in  the  mean  square  sens:)  estimator  of  X  which 


is  a  multiple  of  at  the  parameter  value  X  *  /T.  The 

A 

estimators  X^  are  an  infinite  family  of  admissible  esti¬ 
mators.  It  is  interesting  to  see  that  the  U.M.V.U.  esti¬ 
mator  is  best  only  at  the  parameter  v  alues  X  »  09  or  X  ■  0 , 
Suppose  we  have  prior  information  that  X  $  l. 

Theorem  2.  If  we  have  prior  information  that 
for  the  Poisson  density  that  X  i  l  then 


*3  -l*'*  *  4K 


C7) 


is  an  estimator  with  uniformly  smaller  mean  square  error 

A  At 

that  Xj  and  therefore  X^  is  inadmissible. 

Proof. 


i.s.e.  of  Xj  •  X 


^ 


l2  e  21 


3/2 


n 


n2  J 


<  X  *  Rr>  (X) 

Al 

if  X  $  1. 

4 .  The  Normal  Distribution . 

2 

Case  1.  o  Known. 

Consider  an  estimator  of  the  mean  u  normal  distribution 
with  variance  one.  The  U.M.V.U.  estinator  of  the  mean  is 
x  with  variance  1/n.  The  value  of  k  which  minimizes  the  mean 
square  error  among  all  estimator:;  of  the  form  kUj  is 


1/n  ♦  y 


Therefore  an  estimator  of  the  form 


“  T7'n'l>  T  x  where  1  -  0  C9) 

a 

is  the  best  estimator  which  is  a  multiple  of  y^  at  the 

2 

parameter  values  y  ■  l.  It  is  interesting  to  note  that  the 
U.M.V.U.  estimator  is  best  only  at  the  parameter  value  «*. 
Every  one  of  the  above  estimators  is  admissible. 

Suppose  we  have  information  that  the  mean  is  bounded, 
that  is  y2  <  l. 

Theorem  3.  If  we  have  prior  information  that 

2 

y  <  1  then  2  is  an  estimator  with  uniformly  smaller  mean 
square  error  than  Gj  and  therefore  y^  is  inadmissible. 

Proof . 

m.s.e.  of  y  -  *2/n  ♦  y2/n2 

2  l2  ♦  2Z/n  ♦  1/n2 


■  l/n 

<  1/n 


i2  ♦  y2/n2 
l2  ♦  21/n  ♦  1/n2 


R*  (y)  if  y  <  t. 

yi 


Case  2.  o  Unknown. 

With  the  variance  unknown,  x  *  ^  is  the  U.M.V.U.  esti- 

2 

aator  of  u  with  variance  o  /n.  The  value  of  k  which 


aininizes  the  mean  square  error  among  all  estimators  of  the 
fora  kQj  is 


k  - 


o2/n  ♦  u2  1 


ns 


♦  1 


where  S2  « 

a 


(10) 


Consider  an  estimator  of  \i 


V2  "  T7 Tni)  ♦  1  M1 


(ID 


The  estimator  u2  is  best  among  estimators  of  the  form  k 

2  2 

along  the  lines  la  ■  u  in  the  parameter  space. 

Theorem  4.  Suppose  we  have  prioir  information 


that  S  l  then  u2  is  an  estimator  with  uniformly  smaller 

O  A  A 

■ean  square  error  than  Uj  and  hence  is  inadmissible . 
Proof . 


i.s.e.  of  p. 


I  *2g2  ♦  U2/n 

”  i2  ♦  ii  ♦  l/n: 
n 


a_ 

n 


i2  ♦  4 1 
a1  n 

g 2  +  2_£  +  1. 

" 


if  4  <  t. 

o 

5.1^  The  Bivariate  Normal  Distributi on . 

Case  1.  Covariance  Matrix  I. 

Consider  a  bivariate  normal  with  variance-covariance 
aatrix  I  and  mean  vector  unknown.  The  U.M.V.U.  estimator 
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of  y  «  (y  ,  y  )  is  y  *  (X,Y).  It  has  been  shown  to  be 
x  y  i 

admissible.  The  value  of  k  which  minimizes  the  mean  square 

A 

among  all  estimators  of  the  form  ky^  is 


The  estimator 


(2/n)  ♦  ||y|| 


is  best  in  the  mean  square  sense  among  estimators  of  the  form 
kpj  for  points  on  the  sphere  y^2  ♦  yy2  »  1  in  the  parameter 
space.  Again  the  U.M.V.U.  estimator  appears  best  only  at 
the  points  in  the  neighborhood  of  the  point  at  infinity. 

Theorem  5.  If  J |y| |  <  t,  1  >  0  then  u2  is  an 
estimator  with  uniformly  smaller  mean  square  error  than  yj 
and  hence  y^  is  inadmissible. 

Proof. 


I2  ♦  lt*H 

i.s.e.  of  y  *  - - - 

2  n  2  4£  4 

1  +  TT+  ' 7 

ft 


<  =•  »  R- 


if  llJll2  S  t. 


2  2 

Case  2.  Covariance  matrix  1  a  ,  a  Unknown. 

The  U.M.V.U.  estimator  for  the  mean  is  again 

with  variance  2o2/n.  The  value  of  k  which  minimizes 
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-V-- 


NON-NULL  DISTRIBUTIONS  OF  THE  LIKELIHOOD  RATIO 
CRITERION  FOR  THE  SPECTRAL  MATRIX  OF  A 
GAUSSIAN  MULTIVARIATE  TIME  SERIES 


IC  is  well  known  chut  a  consisecne  estimator  of  a  is  given  by  (see 

J  * 

Parzen  (1966),  BriUinger  (1974) 
l  "l 

<Jlk(id)  -  —  £  I  .  (w  ♦  9  )  • 

J*  n  Jk  m 

a  *-n^ 

wh«r«  n  •  2n^  ♦  1,  btlng  th«  bandwidth  parameter.  This  e9tlaate  a  ^(ui) 
can  be  alternatively  written  as  (see  Priestly,  R^p  and  Tong  (1973)J 

V“>  ■  i  \  V“'  ♦  V  • 


whsra  id'  -  id  -  (2n  (n.+l) ] /T.  The  spectral  matrix  S(w)  is  then 

X  X 


S(u)  -  ((Of  (u*)))  -  -  K  say 

X  Jfc  O'* 


.-.(1.1) 


Goodaaa  (1963)  showed  that  A(u»)  is  distributed  as  the  Complex  Wlshart 

x 


distribution  CW(A  ;  p,  n,  £)  defined  by 


'  «(*;  P.  a.  ^)  -  (rp(n))'1  |  E  (0i)rn  |A(tli)|n_l>  exp  (-tr£  “  A)  ,  ...(1.2) 

■v 

where  T(-)  1.  defined  in  James  (1964).  It  Is  easy  to  see  from  (1.1)  that 

► 

S(<u)  is  distributed  as  CU(S  ,  p,  a.  t  ).  The  study  of  the  structures  of 

tha  above  spectral  matrix  and  especially  testing  various  hypotheses  concerning 

this  matrix  arise  in  the  analysis  of  the  data  in  numerous  areas,  such  as  the 

vibrations  of  the  airframe  structures,  meteorological  forecasts  and  signal 

detection  (see  Hannan  (1970),  Liggett  (1972,  1973),  Priestley,  Rao  and 

Tong  (1973),  BriUinger  (1974)  and  Krishnstah  (1976)). 

In  this  paper,  we  consider  the  problem  of  cesting  the  hypothesis 

R  :  -  r(u)  •  E  ((*i)  ,  where  Z0(td)  is  a  known  macrlx  for  all  w  against 
9  »  V  * 

the  alternative  H, :  E  (<d)  i  Ea(w).  We  shall  obtain  the  hth  moment  of  tha 
L  X  V* 

likelihood  ratio  criterion  and  use  these  moments  to  derive  some  represen¬ 
tations  of  the  non-null  distributions  of  the  test  statistic.  These  are 
easily  programmable  on  a  computer  to  obtain  the  exact  power  as  well  as  the 


percentage  points  to  any  degree  of  accuracy. 


42 


2.  SOME  PRELIMINARIES 


The  following  results  and  definitions  are  needed  In  the  subsequent 
sections.  The  Mellin  integral  fransforn  of  a  function  f(x)  of  a  real 
variable  x.  defined  only  far  x  >  0  Is 

M((f(x)|s)  -  E(x*  *)  -  x*  lf(x)dx  ,  ...(2.1) 

where  s  is  any  coaplex  variate.  Under  suitable  restrictions  (Tltchaarsh  (1948)) 
satisfied  by  all  density  functions  considered  in  this  paper,  there  Is  an  Inverse 
formula  or  lnversa  Mellin  transform 

f (x)  -  (2xl)_1  x**{Mf (x) | s Ids  ...(2.2) 

valid  alaosc  everywhere.  A  path  of  integration  is  any  line  parallel  to  the 
imaginary  axis  and  lying  within  the  strip  of  analyclty  of  M(f(x)|s). 

Leans  1.  Lee 

♦  (O  ■  p(x)dx 

be  the  sonant  function  of  e  random  variable  x  with  the  distribution  law 

P(x).  tf 

♦(e)-0(t'k)  ...(2.3) 

with  real  pare  of  e  tending  to  then  9(e)'  can  be  expanded  aa  a  factorial 
series  of  tha  fora 

*(t)  "  C-n  T ( e+s)/r (t-ek+iyee) ,  ...(2.4) 

Q*u  a 

a  being  any  arbitrary  non-negative  constant  (see  Hair  (1940)  • 

l«as  2.  let  Che  asymptotic  series  a^x^  converge  to  the  function  g(x) 
la  the  neighborhood  of  x  •  0  (or  be  Its  asymptotic  expansion  when  x  *  0) . 

We  then  have 

e,(l0  -  1  ♦  t"-l  6^  ,  . . .  (2.  i) 

where  the  coefficients  Sj  satisfy  the  recurrence  relation 

Vj  ZL  k  “k  Vw  flo  - 1  •••<2-6) 

(see  Kalinin  and  Shalaevskll  (1971). 
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Consider  now  the  test  for  the  null  hypothests  H  :  Z(h>)  »  TJ(u),  where 

0  X  \ 

J0(u)  Is  a  known  matrix  against  tha  alternative  H,  :  Z(ui)  d  Proceeding 

"V  1  A#  'X.0 

aa  la  Anderson  (1958,  p.  264),  It  can  be  shown  that  the  likelihood  ratio 


criterion  la 


X  -  (e/n)pB  | A  r0'1  |“  etr  (-A  E0"l>  . 

where,  for  convenience,  we  have  suppressed  ut  and  written  A  for  A(j),  etc. 

-v,  -v 

Alternatively  X  may  be  expressed  as 
X  -  (e/n)pn  |A|"  etr  (-A). 

where  A  %  CW  (A,  p.  ^  ).  Et  -  z  Z^1  . 

Lean*  J-  Let  z  :  pxp  be  a  complex  symmetric  matrix  whose  real  part  Is 
positive  definite  and  S:  pxp  be  a  hermltlan  matrix.  Then 


r . 


3  _S  >  o 


exp(-tr  Z  S)  |S|  4~p  dS 

■V.  A.  A.  A, 


*  r.  (n)  \1\' 

9  v 


(asm  mg.  9*  ot  James  (1964)), 


The  h  ooaent  ot  the  statistic  X  defined  la  (2.8)  la  given  by 

‘a,  a JT  r  (n) 

P 


Proof:  E(X  )  -  / 


A-  A*>0  CW(*‘  P*  "*  -X  * 

A.  A.  A, 


-  (a/n)pnh  l£ll' 


'  .  W 

A -A*  >  0  ^ 

A.  A.  A, 


o+nh-p 


etr{-(^  ♦*!)£>.  d£ 


Using  Lemma  ),  the  result  fo* lows 


Leva  5.  The  following  Identity  Is  true: 


|I  ^  h  c  | -<»>♦*> 
a  -vl 

•  a  ♦  Vp(h4,n)  |r1|-(h4«)  jr  .  ( i >  'l  a  _  £.-*) | 


•  «*  ;>  ■p<h*">  Co1,  <l-);  vg 


..(2.11) 


(1  +  -)  *  k! 
n 


.  ^ 

where  M  •  I  -  E,~l  and  C  (■)  Is  the  zonal  polynontal  as  defined  In  Janes  (1964). 
\  K 

This  Is  valid  since  h  can  be  chosen  so  thatjl  -  ch  max  tj'^l  <  1  ♦  ~  • 


leva  6.  The  following  expansion  for  the  game  function  holds: 
log  r  (*♦««)  -  log  (2W)**  ♦  (x+h-l/2)log  x  -  X 


-  z  <-»r  »fn(h)  ♦iBfl(«). 


r(r+l)xr 

where  Rb(x)  la  the  reminder  such  that  |  R— (xt)  |  £  _9_  ,  8  being  a  constant 


Independent  of  x  and  Br(h)  Is  che  Bernoulli  polynomial  of  degree  r  and  order 


J.  NOR-NULL  DISTRIBUTION  OP  1  AS  A  MIXTURE  OP 
INCOMPLETE  BETA  FUNCTIONS 


la  this  section,  we  shall  derive  the  non-null  distribution  of  L  •  Xn, 


where  X  Is  given  la  (2.8).  Using  (2.10),  che  hth  moment  of  L  under  the 


alternative  hypothesis  la 
n 


E(Lh)  -  E(Xn) 


-  (e/n)ph  ?  ("*h>  |C  |  h  .  |l  a-  *  E.|-(h-H.) 
ST -  V  \  n  1,11 


..(3.1) 


where  r.«J£  Using  the  Inverse  Mellln  cransfore  (see  (2-2)),  the  density 
'  %  . 


function  of  l  Is  given  by 


^vv.viv^.vv\:v/. 


...  ■ 
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.c+t-  L-h-i  (e/n)Ph  ^  (^,  !Cilh.  II  ♦  j£l|-(h4">  dh  • 

We  now  us#  th#  Identic/  (2.11)  of  Lemma  5  end  write  t  •  n+h  to  obtain 
f(U  -  Ktp.n.Ei.  L) 

Co  ZX  {£  P<L> 

k! 

where 

H  -  I  -  E.  . 

p(L)  -  (2ffl)_1  /c+1~  -L'C  6(0  dt 

c-1" 

♦(t)  -  (#/Opt  e'k  (lp  r(t  +  k  i  -  a) 
a-l  ® 

and 

K(p.  n.  rx.  L)  -  np-A  (f (n+l-a) |Cj_|_B  (e/n)-pn  Ln‘l  . 

B y  leeae  6,  log  6(t)  oay  be  expanded  a# 

log  6(c)  -  log  (2*>p/Z  *  log  t"*  *■  Cl(,r/V) 

where  the  constants  qr  are  given  by 

0f  -  <-Dr‘l  f  Cl  8r+l(Vl-a)  >  r(l*1)} 

and  *  “  pZ/Z^  Thus  we  have 

♦  <t>  -  t""<2*)p/2  (l  ♦  Cl(lr/tr>}  * 

where  the  coefficients  can  be  recursively  computed  using  the  relations 
(2.6)  and  (1.6).  In  ocher  words,  the  coefficients  If  satisfy  the  recurrence 


relation 


..0.9) 


Now  from  (3.7)  and  (3.4)  we  get 

p(L)  -  Out)'1  (2w)p/2  /£ L_t  t"v(U-  C1(ytr)}  dt 

The  Integral  on  Che  right  hand  side  of  (3.9)  can  be  easily  computed  If  v  Is 
an  integer  and  Its  value,  by  Cauchy's  theorem  of  residues  is  2nl  times  the 
sum  of  che  residues  ofLCcv(l+  l^/c*)  at  t  ■  0.  This  Is  easily  seen  to  be 
(7*1)  (2*)p/2  £*-(J  ((-log  L)r*V'1  tr/r(v+r>},  lQ-l 

and  thus  from  (3.9) 

p(U  -  (2n)p/2  £*_0((-log  L)V+C_l  tf/r(v+r))  .  ...(3.10) 

Then  from  (3.3),  che  probability  Chat  L  Is  less  than  any  value  Is 

P(L  i  V  -  yp.  ».  £l)  r".0  Zt  £  CK  (M)  /Q°p(L)dL,  ...(3.11) 

where 

CjCp.  n.  tp  -  (r(nel-a) lEj-0  (e/n)‘pn. 

For  computational  purposes  we  lee  t^n-1  and 

Wl.m  V  -  '<?  *•"  <-!••  t)^-l/r(^)  dL  ...(3.12) 

Integrating  (1.12)  by  pert*.  1C  can  be  easily  shown  that  the  following 
recurrence  relation  holds: 

io„u<v  - 

Wl.mV  ■  (-lo«  V^*l/  r(^r>  *  W -<>•“) 

with  this  notation,  (3.11)  becomes 

pai  V  *  Vp.  id  Co  ek  ii  W- 


L  2 


T 


-  £“.0  Rt  r(t>a)/r(t+v+aM)  ,  ...(3.16) 

where  a  1*  any  arbitrary  positive  constant  and  can  be  chosen  to  govern  the 
rate  of  convergence  of  the  resulting  series.  The  coefficients  Rj  can  be 
explicitly  determined  as  Is  done  below. 

Applying  t erase  6  to  each  of  the  gamma  functions  In  (3.16)  we  get 
log  {r(t+a)/r(t+a+v*D }  •  -<v+l)  log  t  +  (A^/t^) 

where 

At  -  -  BJ+i(w*a+l)| /J(J+1)  .  ...(3.17) 

Thus  we  have 

-(v+1)  ^  • 

r(W-«)/r(c+w+«4.1)  -  Z  U  «■  ,  ...(3.18) 

where  by  Leone  2  and  (3.17),  can  be  recursively  computed.  Using  (3.18)  on 
the  right  hand  side  of  (3.16)  we  get 

♦  Cl  <V*r»  •  Co  e'(,rfl>  V1  +  £j-l(Clj/tJ)  •  ...(3-19) 

equating  the  coefficients  of  like  powers  of  c  on  both  sides  of  (3.19),  It  is 
easy  to  check  the  following  explicit  relations  to  determine  : 

Co  VjCl-y  "  ll  <l-l*2*3-”>  ...(3.20) 

where  lj  *  1  and  *  1. 

How  using  (3.16)  in  (3.9)  and  noting  that  term  by  term  Integration  Is 
valid  since  a  factorial  series  Is  uniformly  convergent  In  a  half-plane 
(see  Doetch  (1971)),  we  have 

p(L)  -  (2w)p/*  E*.0  Rj(2wl)'1  /££  L_t  (T(t+a) /T(t+a+v*j) }  dt  ,  ...(3.21) 

We  now  use  the  well-known  Integral  (Tltchmarch,  1948) 
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(2*1)”  /*“ . J"  x”*  (r(s+a)/r(s+a+v+J) }  ds  •  x  (L-x)  J  t 

C- 1*  ■  .  ■■■  — 

rcv+j) 

and  obcaln 

p(L)  -  (2*)p/2  Zm  R}  L*(l-L)w4'J'1/r(v4-j). 

Then  from  (3.3),  we  have  Che  fallowing  exact  noncentral  c.d.f  of  L  as  a 


■ixture  of  lncoisplete  beta  functions? 


PU  '<  Lq)  -  ^(p.  n.  £1)  Z 


"  j  L 

k-0  K  k! 


■v. 

C 

K 


(M)  . 


0t  Bl  (n*a  ,  v+D/rfv+l),  ...(3. 

where 

Dt  -  B(n+e,  vK) 
and 

•s(P.q)  -  (Xp.q))"1  /*  x>~1  (l-x)q_1  dx 

Is  the  Incomplete  bets  function. 


As  s  particular  ease,  putting  Z  -  ZQ,  l.a.,  M  «  0  ,  we  get  from  (3.22)  the 
null  distribution  of  l  u  the  following  mixture  of  Incomplete  beta  functions: 

»«•<  V  -  v**.  *•  l  >  £t-o  Di  \  (n ...o. 

where  D*  la  the  coefficient  with  k^  -  0. 


4.  ROM-NULL  DISTRIBUTION  OP  1  AS  A  MIXTURE  OP 
INCOMPLETE  GAMMA  FUNCTIONS 

Let  Lj  “  -2q  log  1,  whara  1  la  as  In  eq  (2.8)  and  q  le  an  adjustable 
constant  which  can  be  chosen  to  govern  the  rate  of  convergence  of  the  resulting 
series  and  0  <  q  <  «.  If  C(t)  Is  the  characteristic  function  of  L^,  then 


from  (2.10)  we  have 

C(t)  -  (n/e)”qplt  ?p  (n(l-2qlt) }  (?p(n)}'1 


C(c)--  l^r"  (rp  (n))*'  K()  H(t)  .  ...(4 

k! 

where 

H(t)  -  (c/n)'np<lU  r  (n(l-2qlt).  K»  (1  -2c It )  lt >  .  .  .  .  (« 

P 

where 

rp  (a,K)  -  (e)K  ?p(a) 
ai  defined  In  Constantine  (1966). 

Using  Lemae  6,  we  obtain  after  soae  slnpllf lcatlon 
log  H(t)  *  ^  log  ’f  *  k  i.ig  n  *  pn  logfn/e) 

-  v  log  (n(l-2qlc)}  ♦  $r  (l-2qte)  r  •  ...(' 

where  v  *  p2/2  and 

%  -  <-l)C_1  n’r  [Ej.j  Br+l  (kj  ♦  i  -j)]  (r(r-*-l) )_1  .  ...(< 

to  that 

H(t)  -  <2*)p/2  nk  (n/e)P"  {o(l-2qlt)rV  (1  ♦  ^.^a-Zqlt)"1')  . 

where  the  coefficients  Af  can  be  recursively  computed  using  (4.5)  and  Lemna  2 
Noting  that  (1  -  Sit)*®  Is  the  characteristic  function  of  the  gasae  density, 

ga(B.»)  -R“r(a)}  a®"1  t~*/S,  *  >  0.  ci  >  0.  B>0 

wo  have  fron  (4.2)  and  (4.6)  the  density  of  ■  -2qlogA 
...  .  _  ,  r  a  r"  r  «  ?_  (H) 

E*^  Ar  n‘T  gr+T  (2q,  tx)’  A0  -  l  •••< 

where 

Kz(p.  n.  Zi)  -  (2tr)p/2  (n/e)pn/2  |^|B  np_t  (r(rH-l-a)  )_l  ...( 


and  the  noncentral  c.d.f  of  Is  then  the  following  mixture  of  lncoeplete 
gamna  functions 


t?*?*  WA  VAV"  " 


t»Vi',"',  ■  -'  t. »v  !■■-  r. 


'  \ 


P<4  *  x)  •  K2(p,  n.  £t)  e“_0  Ek  n  CK 


(M) 


kl 


I.  ..  A  n 

r«0  r 


cr4-v  <7q.  *> 


whara 

C  (8.0  -  (8ar(a))"1  r  *  (S.y)  dy 
,  a  x  a 


...<4.«) 


Remarks 

(1)  Taking  £  •  E o  •  or  “  0  la  (4.9),  the  e.d.f  of  L,  «  -2q  log  X  In 

%  A#  ^  \  1 

cha  cantral  casa  Is  tha  following  alxtura  of  Incomplete  gamma  functions 

P(Lt  >x)-K2  (p,  n.  I)  £*  A'  n'V  (2q.  a).  ...(*.10) 

whara  A'  ta  tha  coafflclant  Af  with  kf  «  0. 

(2)  Taking  q»l,  wa  obtain  tha  raprasantatlon  of  tha  e.d.f  of  In  tha 
caatral  eaaa  aa  a  alxtura  of  chl-squara  distributions. 

Tha  raprasantatlon#  of  tha  e.d.fa  of  tha  llkallhood  ratio  criterion 
A  obtained  In  this  paper  aaabla  ona  to  coaputa  tha  exact  power  of  tha  testa 
using  only  tha  tables  of  zonal  polynomials  and  tables  of  either  tha  locoaplata 
bats  or  gsM  functions  and  are  In  S  fora  which  can  ba  aaally  programed, 
further,  tha  Introduction  of  tha  adjustable  constants  to  govern  tha  rata 
of  convergence  and  cha  recurrence  relations  (3.8),  (3.13)  and  (3.20)  are’ 
vell-suitad  to  computations  of  power  as  well  as  percentage  points  for  tests 
of  significance. 
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Probability  of  Destruction  of  a  Point  Target  in  Space 

I.  S.  Przemieniecki 

Air  Foret  Institute  of  Technology 
Wright-Pat terson  Air  Force  Base,  Ohio 


Abstract 

This  paper  presents  an  analytical  method  for  the  calculation  of  the  proba¬ 
bility  of  destruction  of  a  point  target  in  three-dimensional  space  for  a  given 
aiming  offset  and  a  spherically  symmetric  Caussian  distribution  of  the  pro¬ 
bability  density  function.  The  resulting  probability  is  presented  in  terms  of 
two  non-dimensional  parameters  expressed  in  terms  of  the  aiming  offset  r0, 
the  lethal  radius  R,  and  the  variance  a2.  The  resulting  spherical  error  prob¬ 
able  is  compared  with  the  corresponding  two-  and  three-dimensional  cases, 
i.e.  the  circular  and  linear  error  probables.  The  limiting  case  of  zero  aim¬ 
ing  offset  is  also  discussed. 

Introduction 

When  the  dimensions  of  a  target  are  small  compared  with  the  effective 
radius  of  the  weapon  within  which  the  target  can  be  destroyed,  the  target 
itself  can  be  represented  by  a  point  in  the  mathematical  model  of  the 
attack.  For  example,  such  a  model  may  represent  the  case  of  an  anti- 
ballistic  missile  (ABM)  with  a  nuclear  warhead  used  to  destroy  an  incoming 
reentry  vehicle.  Here  the  size  of  the  target  itself  is  minute  in  relation  to 
the  lethal  radius  of  the  ABM  warhead  and  consequently  the  target  can  be 
treated  as  a  point  target.  Due  to  systematic  errors  which  occur  in  any 
weapon  system  the  actual  aiming  point  (mean  point  or  detonation)  will  not 
coincide  with  the  target  and  the  usual  randomness  characteristics  of  the 
weapon  will  produce  the  typical  scatter  around  the  mean  point  of  detona¬ 
tion.  The  systematic  error  can  be  accounted  for  by  specifying  the  aiming 
offset  distance  from  the  target  while  the  randomness  error  can  be 
described  by  the  usual  standard  deviation  in  the  three-dimensional  proba¬ 
bility  density  function. 

This  paper  presents  the  analytical  method  for  the  calculation  of  the  proba¬ 
bility  of  destruction  of  such  a  point  target  in  space  for  a  given  aiming  offset 
and  specified  lethal  radius.  The  probability  density  function  (PDF)  describ¬ 
ing  the  weapon  scatter  in  space  is  assumed  here  to  be  given  by  a  three- 
dimensional  Gaussian  distribution  with  equal  variances.  The  probability  of 
target  destruction  is  shown  to  depend  on  two  non-dimensional  parameters 
expressed  in  terms  of  the  aiming  offset  r0  ,  the  weapon  lethal  radius  R,  and 
the  standard  deviation  o  .  Furthermore,  the  paper  defines  the  spherical 
error  probable  (SEP)  as  a  relationship  between  the  lethal  radius  parameter 
R/  o  and  the  aiming  offset  parameter  r0/<7  and  then  compares  it  with  the 
circular  error  probable  (CEP)  and  the  linear  error  probable  (LEP)  .  The 
analysis  includes  also  a  discussion  of  the  limiting  case  of  zero  aiming  offset. 


Evaluation  oj  the  Probability 


The  analysis  of  the  probability  of  destruction  of  a  point  target  for  two- 
dimensional  cases  has  been  treated  extensively  in  the  past  by  various 
authors  and  tabulated  results  are  available  1  ~a  .  This  paper  provides  a 
natural  extension  of  this  analysis  to  three-dimensional  cases  which  may  be 
of  importance  in  determining  the  effectiveness  of  an  ABM  system.  It  will  be 
assumed  that  the  aiming  offset  point  (i.e.  the  mean  point  of  detonation)  is 
at  a  distance  r0  from  the  target  placed  at  the  origin  of  coordinate  axes  as 
shown  in  Figure  1. 


Mean  point  of 
detonation 


Target 


r/W 
'  1 1 
1 1 

_ m_ 

d0  |  I 
'  I  ‘ 

Nj  ! 


r  sin  f  de 


Figure  1.  Target  in  space 


The  aiming  offset  point  ia  assumed  to  lie,  without  any  loss  of  generality,  on 
the  z-axis  so  that  the  three-dimensional  probability  density  function  based 
on  equal  variances  a2  can  be  expressed  as 


/  v  1  _  -x2  y*  (*~r0)2 

p{z'y’z)=w^exp[^~y--^n 

1  r§  r2  .  rr0cos?l 

where  the  relations  r2-x2+y2^z2  and  z  =rcostp  were  used. 

The  probability  of  destruction  of  the  target  is  calculated  from  the  integral 
R  n  2n 

P-  f  f  J p(x  ,y  ,z)r2sin<pdrd<pd$ 


1  rf  [-(’•-’■o)2]  [  -(t-  4-t*0)2 

= - /  !rexPl — — I — 


rr0cos^ 


r2s\n'pdrd(ftd # 


Subsequent  integration  is  possible  with  the  substitution  of  the  following 
variables: 

p-T /  [vTSo]  (3) 

(  The  symbol  p  denotes  the  quantity  /?/  [V2a]  when  used  as  an  integral 
limit  or  as  a  parameter.) 

Po=ro  /  (*) 


u=p-p0 


v=p+p0 


Substitution  of  the  above  variables  into  Eq.(2)  leads  to 


P- -4 —  /  («  *-p0)exp(-uz)du  -  — = —  f  (v  -p0)exp(-v2)du  (7) 

VrrPo  -Po  V7TPo  p0 

which  can  now  be  integrated.  This  results  in  an  expression  in  terms  of 
exponential  and  error  functions,  given  by 

P=P(R/a.r0/a)  (8) 

=2vkr(',cpr(pt'’”>1" 


WrTpo[<»r/  (p+p„)+er/  (p-p0)]] 


Equation  (8)  is  plotted  in  Figure  2  for  a  series  of  values  of  the  offset  param¬ 
eter  r0/  a  .  When  rg/  a=0  (i.e.  when  p0  =  0)  then  Equation  (8)  reduces  to  a 
special  case  of  zero  offset.  Since  pg  appears  in  the  denominator  of  Equation 
(8).  L'Hopital's  rule  can  be  applied  to  obtain  the  limiting  value  when  r0=0. 
This  leads  to 


P(/?/a,0)=^=-p  e*p(-p8)  +  erf  (p) 
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Figure  2 

Probability  of  destruction  of  a  point  target  in 
space  for  a  given  aiming  offset  r0/  a  and  the 
lethal  radius  R/  a 

Spherical  Error  Fh-obable 

For  the  case  of  zero  aiming  offset  r0=0  the  Spherical  Error  Probable  (SEP) 
is  defined  as  the  radius  R  for  which  50^of  all  detonation  points  in  space  will  fall 
within  the  radius  and  the  remaining  50^will  fall  outside  this  radius.  Thus  if  R 
represents  also  the  lethal  radius  any  detonations  within  the  radius  R  will  destroy 
the  target.  This  concept  can  also  be  extended  to  non-zero  values  of  r0.  This  is 
obtained  simply  by  equating  P  in  Equation  (3)  to  1/2,  so  that 

+V£p0[er/  (p+p0)  +  er/  (p-p0)]]  =  1/  2  (10) 

Equation  (10)  gives  a  relationship  between  R/  a  and  r0/  a  for  the  Spherical  Error 
Probable.  This  relationship  is  plotted  in  Figure  3,  which  shows  clearly  how  the 
lethal  radius  R  must  be  increased  with  an  increase  in  the  offset  parameter  r0/  a 
to  achieve  0.5  (i.e.  50;?)  probability  of  target  destruction.  For  r0/ a>  1,  R^r0, 
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Variation  of  linear,  circular,  and  spherical 
error  probable  with  the  aiming  offset  Tq/o  and 
the  lethal  radius  R /  a  . 


For  completeness,  it  is  interesting  to  compare  similar  relationships  for  the 
circular  and  linear  errors  probable.  Thus  using  results  from  Reference  [1]  it  can 
be  shown  that 

£ff„ftn  =  l/2  (“> 

n*  X 

represents  the  required  relationship  for  the  Circular  Error  Probable  (CEP), 
where 

g,=  l-ezp(-p2)  (12) 

1=V»~  (,3> 
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/>n=  frT-l);  Po(w'1)exP(~Po)  (14) 

with  p  and  po  defined  by  Equations  (3)  and  (4),  except  that  R  and  r0  refer  now  to 
the  radii  in  the  two-dimensional  impact  plane  containing  the  target. 

It  can  be  demonstrated  that  the  corresponding  relationship  for  the  Linear 
Error  Probable  (LEP)  for  a  point  target  in  one-dimensional  space  (range)  is  given 

by 


R  — rc 
'/2cr 


-erf 


-(R+Tq) 

V2  a 


=  1/2 


(15) 


where  R  refers  to  the  lethal  range  and  r0  is  the  linear  aiming  offset.  Both  rela¬ 
tionships  for  the  circular  and  linear  cases  are  also  plotted  in  Figure  3. 
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Abstract 

A  new  variational  principle  is  presented  which  enables  the  use  of  a  direct 
approach  to  obtaining  solutions  to  problems  in  elastodynamics  with  mixed  boundary 
conditions.  The  principle  may  be  viewed  as  a  modification  of  Hamilton's  prin¬ 
ciple,  in  which  the  requirement  that  virtual  displacements  do  not  satisfy  the 
displacement  boundary  principle  to  general  elastodynamics.  A  procedure  is  given 
for  using  the  method  to  determine  the  natural  modes  and  frequencies  of  membrances 
and  plates  having  mixed  boundary  conditions. 


I.  The  Variational  Principle 


In  the  application  of  Hamilton's  principle,  it  is  necessary  to  seek  an 
extremal  value  of  the  Lagrangian  function  over  a  class  of  admissible  virtual 
displacements  which  vanish,  for  all  time,  over  the  portion  of  the  boundary 
where  the  displacements  are  to  be  prescribed^".  For  many  problems,  finding  such 
functions  is  not  readily  done;  therefore  it  is  of  interest  to  investigate  the 
possibility  of  extending  the  class  of  admissible  functions  to  include  functions 
which  do  not  satisfy  the  boundary  conditions  on  those  segments  of  the  boundary 
where  tractions  are  prescribed,  nor  on  those  portions  where  displacements  are 
prescribed . 

This  relaxation  of  requirements  on  the  trial  functions  is  precisely  that 

2 

permitted  in  the  use  of  Reissner's  principle  ,  which  is  applicable  to  both  the 

3  4 

static  response  the  the  simple  harmonic  motion  ’  of  elastic  systems.  It  is 
of  interest  to  begin  as  with  Hamilton's  principle  and  to  derive  a  somewhat 
more  general  result. 

Ve  begin  by  writing  a  Lagrangian  function 

L  =  U  -  K  «-  A  (1) 

where  the  atraln  energy,  U,  of  an  elastic  material  is  written  as  a  func¬ 
tion  of  strains, 

U  -=  |  W(e^ j )dV  (2) 

V 

the  potential  energy  of  conservative  external  forces  is  written  as 

A  *  -  {  T1*u1«is  -  |  F^dV  (3) 

S  V 

o 
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and  che  kinetic  energy  Is  expressed  as 


t  .  j  i  p  Vl4V 


The  volume  of  Interest,  V,  is  enclosed  by  the  surface  S  ■  Su  ♦  So»  where 

S  is  the  portion  of  the  boundary  on  which  displacements  are  to  be  pre- 

* 

scribed,  and  S  is  the  portion  on  which  tractions,  T,  ,  are  given.  The 
body  force,  Fj,  is  presumed  to  be  prescribed.  The  customary  requirement 
that  the  displacements  satisfy  a  boundary  condition 


u.  m  u,  on  S 
i  i  u 


is  replaced  by  a  constraint. 


\  »  |  Vui  "  u4*) 


where  the  are  Lagrange  multipliers.  We  assume,  as  in  the  further 

5 

generalisation  of  Relssner's  principle  due  to  Vashlzu  ,  that  the  strains 
and  displacements  may  be  varied  independently.  Thus,  the  requirement  of 
satisfaction  of  the  strain  displacement  equations  is  replaced  by  a  con¬ 


straint 


C2  *  |  ‘ijOij  -  I  (“l,j  *  °J.i,:1,,V 


where  the  comma  denotes  partial  differentiation.  The  strains  are  assumed 
to  be  symmetric;  thus  the  are  also. 

The  new  functional  is 

L*  -  V  -  K  ♦  A  -  Cx  -  C2  (7) 

and  may  be  recognized  as  a  time  dependent  version  of  the  functional  used 
In  Washizu's  generalization.  We  now  seek  to  determine  conditions  under  which 
the  time  Integral  of  the  modified  Lagranglan  function  assumes  a  stationary 
value,  or 
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f 2  *,  . 

J  L  <V  V  *ij-  V  ri 


)dt  =  0 


with  all  21  arguments  of  the  integrand  varied  independently.  If  the  trial 
functions,  ,  and  the  X  have  sufficient  continuity  as  to  permit  the 
necessary  application  of  the  divergence  theorem,  and  if 

•iijUj)  =  «u1<t2)  =  0  (9) 

we  find  the  vanishing  of  the  first  variation  necessitates  that 


I  i'tb-  Ws'ijdvdc  *  1 2 1  *-» 

*  V  #■  M 


*  pu^  -  F^Jfiu^dvdt 


C1  V 


h  V 


2 

- \  \  'eij  -  \  <ui,j  ♦  u1>i),Mtidy,,t  *  1  I  ‘Yu  -  rt)‘Vs,,t 


ci  v 


t.  S 
1  u 


I  I  *VJXiJ  ~  Ti  *<uids<lt  “}  )  ^ui  ~  u*}6T1dsdt  =  0 


t,  S 
1  a 


t.  S 
1  u 


We  recognise  from  the  first  integral  that  the  Lagrange  multipliers 
X^j  have  the  physical  interpretation  of  the  components  of  stress  in  an 
elastic  body,  °£j»  *nd  from  the  last  that  the  have  the  physical  inter¬ 
pretation  of  the  components  of  traction,  T^,  since  the  Vj  are  the  com¬ 
ponents  of  the  normal  vector  at  a  point  on  the  surface.  Thus,  three  sets 
of  Euler  equations  result 


3V  4  xr 

- -  *  a,  ,  in  V 

(11) 

*®ij  ij 

•• 

u.j  *  ri  •  *"i  *"  v 

(12) 

*  I  (“l,J  *  "),!>  In  V 

(13) 

and  the  necessary  boundary  conditions  are  seen  to  be 
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i  J  ij  U  0 

I  ' 

T.  -  T  *  on  S  ( 

i  i  o 

u  ■  u  *  on  S  ( 

1  i  u 

Since  these  are  the  field  equations  and  boundary  conditions  of 
elastodynamics,  we  conclude  that  the  proposed  extension  is  appropriate, 
whether  it  be  viewed  as  a  relaxation  of  the  class  of  functions  to  be  used 
with  Hamilton's  principle,  or  as  an  extension  of  Relssner's  principle  to 
other  than  simple  harmonic  motions.  Equation  10,  it  should  be  noted, 
is  included  in  the  further  generalization  due  to  Yu^  .  This  dynamic 
variational  principle  differs  from  that  of  Dean  and  Plass  7  in  that 
6u  need  not  vanish  on  S..  - 


II.  An  Approximate  Method 

Our  interest  here,  of  course,  is  not  in  deriving  once  again  the 
equations  of  elastodynamics ,  but  rather  in  developing  a  procedure  whereby 
approximate  solutions  can  be  obtained  for  otherwise  Intractable  problems. 

To  do  this,  we  return  to  equation  10  and  make  use  of  the  identification 
that  the  Lagrange  multipliers  are,  in  fact,  stresses  and  tractions. 

We  are  Interested  in  the  class  of  problems  for  which  a  number  of 
solutions  to  the  field  equations  are  readily  obtained,  but  for  which  dif¬ 
ficulties  in  obtaining  solutions  arise  because  of  mixed  boundary  conditions 


in  the  form  of  Equations  15  and  16.  We  assume  that  a  large  number,  N,  of 
systems  of  stresses,  c^n,  *n<*  displacements  u^n,  can  be  found,  and  that 
each  system  satisfies  Equation  12,  and  further,  that  these  stresses  sat¬ 
isfy  an  elastic  constitutive  relationship,  Equation  11,  together  with 
strains  e^”  which  are  derived,  through  Equation  13,  from  the  displacements 
Uin.  We  then  take  as  trial  functions  the  superpositions  of  such  solutions. 


V*'0  '  J, 


and 


N 

u.(x,t)  *  l  anu1n(x,t)  (18) 

n=l 

It  follows  that  every  such  trial  function  will  also  satisfy  the  system 
field  equations,  11-13.  Tractions,  T^,  may  be  constructed  from  these 
stresses  so  as  to  satisfy  Equation  14  on  each  point  of  S. 

For  such  trial  functions,  all  but  the  last  two  Integrals  of  Equation 


10  are  sero  identically.  The  variational  principle  has  thus  led  to  the 
requirement  that 


(19) 


S 

A  as, 


#3 


|  I  (T i  ~  Ti*)  auidsdt  |  (ui  "  ui*)«T1dsdt  -  0 


t,  S 
1  a 


t.  S 
1  u 


From  this  condition,  we  propose  to  construct  en  algorithm  for  determining 

the  set  of  coefficients  *n  which,  for  any  N  selected,  leads  to  the  best 

approximate  solutions  in  the  form  of  equations  17  and  18.  The  most  general 

arbitrary  variation  within  the  space  spanned  by  these  N  solutions  may  also 

be  written  as  an  expansion  of  these  sane  N  solutions,  or 

N 

do  «  l  da  o  “(x,t)  (20) 


as 

du  =  l  da  u  ®(x,t)  (21) 

1  «  ID  X 

m  si 

Substituting  these  and  Equations  17  and  18  into  equation  19  leads  to  a 
single  equation.  Since  the  variation  must  be  arbitrary,  however,  any  con¬ 
venient  choice  of  coefficients,  6*a>  may  be  made.  It  is  particularly  con¬ 
venient  to  make  the  selection: 

da  «  0  for  m  4  p,  m  *  1,N  (22) 

18 

da  *  1  for  m  =  p,  p  *  1,N  (23) 

ID 

so  as  to  obtain  a  system  of  N  equations  for  the  N  coefficients  aR,  or 

i  •» l  (l  V)iVd* - 1  -.V,,'"* 

*1  *.  S 

1  °  u  (24) 

t2 

■  J  ‘I  TiVd*  -  J  “t*  YjiP<,*,dc 

t,  s  s 

la  u 

It 

For  any  given  boundary  conditions,  ^  and  ut  ,  and  for  any  selected  set 
of  solutions  to  the  field  equations,  the  a^  may  be  found  and  the  approximate 

g 

solution  determined.  It  has  been  previously  found  that  static  problems 
with  mixed  boundary  conditions  may  be  successfully  treated  in  a  similar 


manner. 


The  procedure  may  be  applied  to  construct  an  approximate  steady  state 
response  to  a  prescribed  boundary  excitation 


★  ★ 

=  6  (s)cos  Qt  on  S 

1  o 

(25) 

(26) 

=  U.  (s)cos  Ct  on  S 

1  u 

We  require  N  solutions  to  the  field  equations,  of  the  form 


The  time  Integrals  have  been  eliminated  from  Equation  26,  since  the  equa¬ 
tions  must  be  satisfied  at  every  Instant.  Postmultiplying  by  the  Inverse 
of  the  array,  C,  completes  the  algorithm  for  computing  the  coefficients, 
an«  For  Su  *  0,  this  procedure  reduces  to  that  previously  used  to  deter¬ 
mine  the  response  of  an  elastic  strip  to  a  time  harmonic  end  loading*®. 

In  order  to  construct  an  approximate  solution  for  the  free  vibration 

* 

of  an  elastic  solid  with  homogeneous  boundary  conditions,  T^  =  0  on 

* 

and  *  0  on  Sy  We  again  require  that  N  solutions 


■  S^jn(w,x)cOS  wt 


U.”  *  U.n(«,x)co*  wt 
1  1 


be  available.  Here,  Che  natural  frequency,  w,  is  Co  be  determined.  Equa¬ 
tion  24  now  reduces  to  a  set  of  homogeneous  algebraic  equations 


7  a  C  (u)  =  0  (34) 

n=l  "  "P 

with  C  given  by  Equation  30.  Estimates  of  natural  frequencies,  u,  arise 
np 


from  setting 


det(C  )  =  0 
np 


and  the  corresponding  mode  shapes  may  be  determined  by  returning  the  re¬ 
sulting  value  of  o)  to  Equation  34  and  computing  ratios  of  coefficients. 

Ve  have  new  completed  the  development  of  a  procedure  whereby  approx¬ 
imations  to  the  forced  response  of  elastic  solids  to  time  harmonic  boundary 
excitations  may  be  obtained.  Approximations  to  the  natural  frequencies 
and  the  corresponding  mode  shapes  may  also  be  developed,  even  for  objects 
of  complex  geometry.  If  the  available  set  of  solutions  can  be  shown  to 
be  complete,  convergence  to  an  exact  answer  is  to  be  expected  in  both  cases. 
If  the  available  set  of  solutions  Is  merely  "large",  only  an  approxlmar ion 
can  be  anticipated. 
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III.  Summary 


A  new  variational  principle  has  been  developed  and  used  to  develop  a 

general  approach  to  the  problem  of  determining  the  frequencies  and  mode  shapes 

of  freely  vibrating  elastic  systems.  Specific  results  from  an  application  to 

9  10 

the  vibration  of  membranes  have  been  given  elsewhere  ’  .  These  results  may, 

of  course,  be  immediately  transferred  to  other  physical  systems  having  the  same 
differential  equation  and  boundary  conditions,  such  as  the  waveguide.  More 
recently,  the  method  has  been  applied  to  the  determination  of  the  mode  shapes 
and  natural  frequencies  of  vibrating  plates.  Specifically,  the  problem  con¬ 
sidered  was  the  circular  plate,  clamped  on  a  portion  of  the  edge  and  free  on 
the  remainder*^.  Even  though  the  solution  algorithm  requires  finding,  by 
iteration,  the  values  of  a  parameter  which  zero  the  determinant  of  a  large 
matrix,  the  computation  time  required  was  found  to  be  substantially  less  than 
that  required  to  obtain  similar  results  by  using  the  method  of  finite  elements. 
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ABSTRACT 

The  linear  velocity  and  velocity  profile  of  argon  flowing 
in  a  flow  tube  were  determined  experimentally  by  use  of  a 
luninosity  time-of-f light  technique.  The  experiment  can  serve 
as  a  simple  introduction  to  flow  tubes  and  to  velocity  profile 
determinations  in  a  physics  or  fluid  mechanics  laboratory.  The 
experiment  allowed  the  observation  of  a  laminar  flow  and  a  com¬ 
parison  of  linear  velocity  calculated  from  plug  flow  considera¬ 
tions  with  the  linear,  centerline  velocity  determined  experimen¬ 
tally.  It  was  found  that  the  linear  velocity  was  1.62  times  the 
plug  flow  (average)  velocity  at  a  flow  rate  of  4576.0  standard 
cc/min  and  1.77  times  the  plug  flow  velocity  at  5262.4  standard 
cc/min.  The  velocity  results  were  compared  with  those  predicted 
by  the  Poiseuille  flow  formula  which  governs  ideal,  fully  developed 
internal  flows. 

INTRODUCTION 

Flow  tubes  are  relatively  simple  devices  vrfiich  can  be  used 
to  determine  the  rates  of  fast  reactions,  to  generate  excited 
chemical  species,  and  to  study  the  rates  of  decay  of  the  emissions 
from  excited  species.  These  studies  can  be  accomplished  because 
in  the  flowing  gas,  little  mixing  occurs  between  reactants  enter¬ 
ing  and  products  leaving  the  tube.  The  determinations  of  the 
mechanism  and  kinetics  of  a  large  number  of  fast  reactions  at 


temperatures  below  1000  K  have  been  done  in  this  manner1  ’  J .  In 
addition,  with  improvements  in  flow  tube  systems ,  high  temperature 
reactions  occurring  in  such  systems  as  Fe/C^^  and  Al/0^^  have 
been  studied.  Other  applications  of  flow  tubes  have  been  report- 
ad[5'  «. 

In  order  to  perform  meaningful  experiments,  the  linear  velo¬ 
cities  of  the  flowing  gases  must  be  determined  accurately.  Linear 
velocities  have  been  reported  in  many  of  the  published  studies^’  * 
These  velocities  were  needed  to  calculate  reaction  rates  and  to 
provide  a  basis  for  a  correlation  between  various  physical  para¬ 
meters.  In  general,  the  velocities  were  calculated  using  the  plug 
flow  assumption^  .  Plug  flow  is  defined  as  an  idealized  state 
of  flow  such  that,  over  any  cross-section  perpendicular  to  the 
fluid  motion,  the  mass  flow  rate  and  the  fluid  properties,  P,  T, 
and  p,  remain  constant.  Extensive  information  is  lacking  concern¬ 
ing  the  validity  of  the  plug  flow  assumption  in  flow  tube  applica¬ 
tions  where  the  pressure  is  low  (4  -  10  Torr) ,  tut  where  the  flow 
is  still  in  the  viscous  flow  regime.  Swearengen^11^  and  Armstrong 
and  Davis1  1  have  reported  limited  results  on  this  problem, 
although  it  was  not  a  main  concern  of  their  studies.  The  current 
work  was  conducted  to  develop  a  simple  experimental  method  to  de¬ 
tail  the  validity  of  the  plug  flow  assunption  in  the  calculation 
of  the  linear  velocity  of  a  low  pressure  gas  in  a  flow  tube.  The 
linear  velocity  was  determined  experimentally  and  compared  to  cal¬ 
culations  based  on  the  plug  flow  assunption.  In  addition,  the 
velocity  profile  was  studied  to  determine  possible  deviations  from 
plug  flow.  The  experimental  apparatus  vhich  was  utilized  and  the 
technique  for  measuring  the  velocity  profile  are  simple  enough  so 
that  they  could  be  used  as  a  laboratory  experiment  in  a  physics 
or  fluid  flow  laboratory.  The  experiment  cculd  introduce  the  stu¬ 
dents  to  the  use  of  flow  tubes  and  also  could  serve  as  a  simple 
method  of  obtaining  the  velocity  profile  for  a  gas  in  laminar 
flow.  The  equipment  is  inexpensive  and  simple  to  operate,  yet  the 
results  give  a  good  approximation  of  the  theoretical  parabolic 


FIGURE  I:  Schematic  diagram  of  the  flow  tube  system.  The  component  parts  are:  a,  needle  value;  b,  flow 
meter;  c,  ball  value;  d,  trigger  module;  e,  square  wave  generator;  f,  copper  electrodes;  g, 
oscilloscope;  h,  optical  bench;  i,  photomultiplier  tube  (PMT)  system;  j,  x-y  recorder;  k, 
signal  averager;  1,  FMT  power  supply;  m,  anplifier;  n,  vacuum  gauges;  o,  capacitance  manometer 
pressure  gauge;  p,  flow  tube. 


velocity  profile  to  be  expected  under  the  experimental  conditions. 

While  standard  techniques,  such  as  pitot  tubes  and  hotwire 
anemometers^’  ^  have  been  used  in  the  past  for  velocity  pro¬ 
file  measurements ,  the  present  experiment  utilized  a  luminosity 
time-of-f light  technique. 

A  stream  of  argon  gas  was  excited  in  a  high  voltage  discharge 
in  the  flow  tube,  near  the  entrance.  The  time-of-f light  of  the 
luminous  argon  was  determined  at  a  series  of  stations  downstream 
of  the  discharge  by  detection  of  the  visible  emission  of  the 
excited  argon.  The  linear  velocity  of  the  flowing  gas  was  calcu¬ 
lated  from  these  data.  The  advantage  of  this  technique  for  a 
laboratory  environment  is  the  simplicity  of  the  approach.  Once 
the  apparatus  is  assembled,  simple  movement  of  the  detector  along 
an  optical  bench  is  the  only  adjustment  which  needs  to  be  made. 

By  traversing  the  detector  perpendicular  to  the  direction  of  gas 
flow  at  various  stations  downstream  at  the  flow  tube,  the  develop¬ 
ment  of  the  profile  can  be  observed. 

EXPERIMENTAL  SECTION 

The  flow  tube  arrangement  utilized  in  this  experiment  was 
based  on  the  system  described  by  Swearengen,  Davis  and 
Niemczyk^7,  ^ .  It  is  shown  schematically  in  Figure  1.  The  flow 
tube  was  composed  of  sections  of  thick-walled  pyrex  glass  tubing 
5.1  cm  I.D.  joined  together  to  form  a  system  162  cm  in  length. 

The  tube  consisted  of  a  gas  injection  region,  a  discharge  region 
and  a  test  region.  It  was  supported  on  an  aluminum  stand. 

A  gas  flow  system  furnished  the  flow  tube  with  a  continuous 
supply  of  high  purity  argon  ( 99 . 98%) .  Since  the  volumetric  flow 
rate  was  required  in  order  to  calculate  the  theoretical  linear 
velocity  of  the  gas,  this  flow  rate  was  measured  using  a  Hastings 
Model  A11-10K  Linear  Mass  Flowmeter.  The  flow  rate  was  controlled 
by  a  Jamesbury  ball  valve  placed  downstream  from  the  flowmeter. 

The  gas  entered  the  tube  in  the  injection  region.  This  31  cm 
length  of  the  flow  tube,  located  directly  upstream  of  the  discharge 


region,  was  provided  in  order  to  obtain  a  reattached  flow  of  gas 
prior  to  entry  into  the  discharge  region.  Provision  for  the  gas 
flow  to  become  reattached  in  the  flow  tube  becomes  important  when 
the  gas  feed  line  diameter  is  different  from  the  flow  tube  dia¬ 
meter.  That  is,  vrtien  gas  is  suddenly  expanded  as  it  enters  a  tube, 
a  jet  stream  develops.  Therefore,  it  is  necessary  to  find  a  dis¬ 
tance  downstream  from  the  inlet  hole  where  the  flow  becomes  re¬ 
attached.  In  the  present  study,  the  flow  tube  diameter  was  5.1  cm 
and  the  inlet  hole  diameter  was  1.3  cm.  For  this  configuration, 
the  recirculation  zone  length  was  found  to  be  15  cm^5^ .  A  section 
of  tube  approximately  twice  this  length  was  used  to  serve  as  the 
injection  region. 

Once  reattached  flow  was  achieved,  the  gas  entered  the  dis¬ 
charge  region.  A  portion  of  the  argon  flowing  between  copper 
electrodes  was  ionized  and  excited  into  metastable  electronic 
states.  These  states  continued  to  emit  detectable  visible  radia¬ 
tion  for  a  considerable  distance  downstream,  and  this  emission 
(or  afterglow)  served  as  the  basis  for  the  velocity  determinations. 

The  discharge  region  was  constructed  with  a  pyrex  glass 
cross,  20  cm  in  length.  The  cross  was  located  directly  downstream 
of  the  gas  injection  region.  Copper  electrodes,  0.32  cm  in  dia¬ 
meter,  were  inserted  in  the  cross  arms  through  1.3  cm  thich  plexi¬ 
glass  plates.  The  electrodes  were  adjusted  to  provide  a  1.0  cm 
gap  perpendicular  to  the  flow  tube  centerline.  A  high  voltage  arc 
was  produced  between  the  electrodes  with  an  EG&G  Model  TM-11  30  kv 
trigger  module.  A  repetitive,  pulsed  discharge  was  achieved  with 
a  2  Hz,  20  v  peak-to-peak  square  wave  from  a  Wavetek  Model  IV 
voltage  controlled  generator. 

The  observation  section  of  the  test  region  of  the  flow  tube, 

61  cm  in  length,  was  placed  directly  downstream  of  the  discharge 
region.  An  optical  bench  was  positioned  alongside  this  portion  of 
the  tube.  The  bench  contained  a  scale  graduated  in  cm  to  measure 
distances  from  the  electrode  to  the  observation  point.  A  IP21 
photarultipiier  tube  (FMT),  used  to  detect  the  radiation  from  the 


argon  afterglow,  was  mounted  on  the  optical  bench.  A  variable 
slit,  adjusted  to  a  width  of  0.2  cm,  was  placed  in  front  of  the 
PMT.  An  iris,  adjusted  to  give  a  diameter  of  1  cm,  was  placed  in 
front  of  the  slit.  The  iris  allowed  the  photomultiplier  tube  to 
sanple  the  radiation  at  different  heights  in  the  flow. 

The  output  of  the  PMT  was  fed  into  a  Keithly  Model  427  cur¬ 
rent  amplifier.  The  signal  was  then  sent  into  a  Princeton  Applied 
Research  Model  TDH-9  Waveform  Eductor  (signal  averager).  Approxi¬ 
mately  100  signals  were  averaged.  The  signal  averager  was  trig¬ 
gered  by  the  signal  from  the  square  wave  generator..  Therefore, 
the  signal  averager  accepted  signals  only  when  the  discharge  was 


triggered.  The  averaged  signal  was  then  sent  into  a  Hewlett- 
Packard  Model  7045a  X-Y  Recorder.  The  recorder  produced  a  curve 
of  relative  intensity  vs  time. 

A  20  cm  long  pyrex  glass  tee  was  used  in  the  test  section 
following  the  observation  area.  This  tee  housed  the  pressure 
and  vacuun  gauges.  Hastings  DV-4D  and  DV-6  vacuun  gauges  were 
used  to  measure  vacuum.  Pressure  during  gas  flow  measured  with 
a  capacitance  manometer  device  (MKS  Baratron  Type  77  Pressure 
Meter) . 

Two  Welsh  17.7  CFM  mechanical  vacuun  pumps  were  used  to  re¬ 
duce  the  pressure  in  the  flow  tube.  A  pressure  of  50  microns  was 
achieved  with  no  gas  flow.  A  ball  valve,  used  to  retard  the  gas 
flow  to  provide  a  higher  gas  pressure  with  no  change  in  flow  rate, 
was  placed  in  the  line  between  the  pimps  and  the  flow  tube. 

RESULTS 

A.  Linear  Velocity 

The  linear  velocity  of  the  flowing  argon  gas  was  experimen¬ 
tally  determined  at  various  gas  pressures  and  flow  rates  by  de¬ 
tecting  the  afterglow  from  the  excited  argon  with  the  PMT.  For 
these  measurements,  the  PMT  was  positioned  at  various  stations 
dowistream  of  the  electrodes.  The  height  of  the  PMT  system  was 
adjusted  such  that  a  signal  was  detected  from  the  flow  tube  center 


line.  The  pressures,  flow  rates  and  downstream  station  distances 
for  vrtiich  linear  velocities  were  calculated  are  listed  in  Table  1. 
The  centerline  signal  from  the  flow  tube  was  determined  by  travers¬ 
ing  the  FMT  and  the  iris  vertically  at  each  station  until  the  fast¬ 
est  time-of-f light  of  the  glowing  gas  was  obtained. 

Data  gathered  during  this  operation  included  the  distance  (X) 
of  the  FMT  station  downstream  of  the  electrodes  and  the  time  (t) 
at  which  the  FMT  detected  the  leading  edge  of  the  emission  of  the 
argon  afterglow.  An  adjustment  was  made  in  the  value  of  X  in  order 
to  correct  for  the  self-diffusion  of  the  glowing  argon.  The  square 
root  of  the  mean  square  distance  diffused,  (X2)  ,  was  calculated 
using  Eq  (1)[16’  17] 

TX1}*  =  /2Dt  (1) 

where  D  is  the  diffusion  coefficient  of  argon  and  t  is  the  time- 
of-flight  of  the  gas  to  each  station.  The  correction  amounted  to 
less  than  4%  of  the  total  distance  traversed. 

The  linear  velocity  of  argon  was  calculated  by  performing  a 
linear  regression  analysis  on  the  data.  A  linear  least  squares 
fit  was  computed  by  regressing  t  on  X.  Plots  of  t  vs  X  were  con¬ 
structed  along  with  the  lines  representing  the  results  of  the 
least  squares  analyses.  This  information  is  shown  in  Figure  2. 

The  linear  velocities  of  the  gas  flows  were  determined  by  invert¬ 
ing  the  slopes  of  the  lines.  The  results  are  shown  in  Table  II, 
along  with  the  plug  flow  velocity  calculations. 

Upf  =  1.43Q  P(ATM)/60  P(GAS)  A  cm/sec  (2) 

vhere  Q  is  the  volunetric  flow  rate  in  standard  cc/min  of  air,  1.43 
is  the  conversion  factor  for  the  Hastings  Mass  Flowmeter  to  convert 
Q  into  the  volumetric  flow  rate  of  argon,  P(ATM)  is  the  atomspheric 
pressure  in  torr,  P(GAS)  is  the  flowing  gas  pressure  in  torr,  A  is 
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CLn,e~?f_fL?-ShJ  (t)  of  the  glowing  argon  gas.  The  experimental 
conditions  are  given  m  Table  I.  The  data  points  are  designated  by  the 
symbols  O,  0,  A,  and  +  for  data  taken  at  5.0,  4.2,  6.0,  and  10. 1  torr. 
respectively. 


the  cross-sectional  area  of  the  flow  tube  (20.27  cm2)  and  60  is 
the  conversion  factor  from  minutes  to  seconds. 

B.  Velocity  Profile 

The  velocity  profile  of  the  gas  flow  is  dependent  upon  the 
type  of  flow  present  (i.e.,  molecular,  viscous,  or  intermediate). 
The  division  between  these  flow  regimes  is  determined  by  the  value 
of  the  Knudsen  lumber  L  .In  this  study,  the  mean  free  path  of 
argon  at  5.0  torr  was  calculated  to  be  1.04  x  10~*  cm  and  the 
Knudsen  number  was  2.05  x  10  4.  Therefore,  the  flow  is  in  the 
viscous  flow  region1  J .  Furthermore ,  the  Reynolds  lumber  was 
calculated  to  be  280,  which  showed  that  the  flow  was  laminar^ 10  ^ 
as  well.  Therefore,  a  parabolic  velocity  profile  was  predicted 
to  form  in  the  flow  tube. 

The  instantaneous  linear  velocity  of  argon  was  determined  at 
various  heights  above  and  below  the  flow  tube  centerline  from  time 
vs  distance  data.  The  data  for  these  calculations  were  obtained 
at  stations  28.5  cm  and  35.9  cm  downstream  from  the  electrodes. 
These  experimental  velocities  were  then  plotted  against  distance 
above  and  below  the  flow  tube  centerline  to  obtain  the  profile. 

The  results  are  shown  graphically  in  Figure  3.  The  solid  para¬ 
bolic  curve  in  Figure  3  is  a  theoretical  velocity  profile  charac¬ 
teristic  of  a  fully  developed  internal  flow.  This  profile  was 
determined  using  the  Poiseuille  flow  formula  for  a  flow  in  circu¬ 
lar  cross-section . 


u(y)  =  2uav  [1  -  y2/R2) ] 


(3) 


where  y  is  the  radial  distance  from  the  centerline  and  R  is  the 
radius  of  the  tube,  u^  is  defined  in  Eq  (4) 


-R2  &p 
UAV  _  8n  Ax 


(4) 


where  n  is  the  viscosity  and  Ap  is  the  pressure  difference  due  to 


viscous  forces  over  a  length  of  tube  ax.  If  the  pressure  at  one 

end  of  the  tube  is  known,  then  the  pressure  at  the  other  end  can 

r  o  i 

be  calculated  from  Eq  (5)L  J 

P2*  -  Pj/  =  16F'  Ln  RqT/wR'  (5) 

vhere  F"  is  the  flow  rate  in  moles/sec,  L  is  the  length  of  the 
tube,  Rq  is  the  universal  gas  constant,  and  T  is  the  absolute 
temperature.  For  example,  using  for  the  viscosity  of  argon 
223.7  x  10  *  poise  at  296°  K,  a  P^  of  10.1  torr  and  a  flow  rate 
of  3.62  x  10  *  mole/sec,  the  pressure  difference  due  to  viscous 
forces  was  calculated  to  be  1.09  x  10  1  torr. 

The  experimental  velocity  profiles  in  Figure  3  revealed  an 
approximately  parabolic  configuration  at  both  stations.  While 
the  profile  obtained  at  35.9  cm  downstream  slightly  more  closely 
resembled  the  theoretical  curve  than  the  profile  at  28.5  cm,  the 
difference  in  the  profiles  is  within  experimental  scatter.  There¬ 
fore,  the  flow  was  considered  fully  developed  at  both  downstream 
stations.  The  characteristic  length  for  an  internal  flow  to  be- 

r  9i 

come  fully  developed  was  calculated  using  Eq  (6)  1 

\  =  .03  (Re)  D  (6) 

where  is  the  characteristic  length,  Rg  is  the  Reynolds  number 
and  D  is  the  tube  diameter.  The  laminar  development  length  was 
calculated  as  42.7  cm  from  the  point  where  the  flow  became  reat¬ 
tached  at  the  inlet  region.  Thus,  the  flow  was  fully  developed 
at  16.7  cm  downstream  from  the  electrodes.  This  sinple  calcula¬ 
tion  was  strongly  supported  by  a  numerical  calculation  for  the 

[19] 

same  system  vhich  was  performed  by  Rapagnani  and  Davis 

DISCUSSION 

The  velocity  profiles  of  the  flowing  argon  gas  at  various  low 
pressures  were  parabolic.  This  result  indicates  that  a  fully 


developed  laminar  flow  was  present  throughout  the  experiment. 

Based  on  this  conclusion,  a  comparison  can  be  made  between  the 
experimentally  determined  values  for  plug  flow  (Eq  (2))  and  the 
average  velocities  calculated  by  the  use  of  Eqs  (4)  and  (5). 

These  values  are  reported  in  Table  II.  The  correspondence  between 
values  shows  that  the  plug  flow  velocity  can  be  used  as  the  average 
velocity  in  calculations  for  flow  tubes  when  low  pressure  laminar 
flow  is  occurring. 

The  results  of  the  linear  velocity  calculation  along  the 
flow  tube  centerline  should  be  twice  the  average  or  plug  flow 
velocity  in  the  ideal  case.  This  is  evident  vtaen  y  is  set  equal 
to  0  in  Eq  (3).  A  comparison  between  the  experimental  linear 
velocity  results  at  y  =  0,  (ug),  and  the  plug  flow  values  shows 
that  the  factor  of  2  is  being  approached,  but  experimentally  the 
linear  velocity  is  somevhat  slower  than  predicted  by  theory.  The 
discrepancy  can  be  the  result  of  a  number  of  factors.  Plug  flow 
calculations  do  not  take  into  account  frictional  effects,  since 
this  calculation  is  based  on  an  ideal  flow.  In  this  study,  the 
flow  experienced  frictional  effects  due  to  shearing  stresses  and 
pressure  losses  which  caused  the  deviation  from  plug  flow.  Also, 
in  the  present  case  the  frictional  effects  were  increased  when  the 
gas  flowed  past  the  electrodes. 

Based  on  these  studies,  it  is  concluded  that  a  calculation  of 
centerline  velocity  based  on  a  laminar  velocity  profile  gives  a 
much  better  value  than  a  calculation  based  on  plug  flow.  If 
perturbations  to  the  flow  are  slight,  as  in  the  present  case,  a 
correction  factor  of  1.6  -  1.7  applied  to  the  plug  flow  value  will 
give  a  good  estimate  of  the  centerline  flow  velocity. 
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